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SIMPLIFIED  MODEL  FOR  PREDICTION 
OF  NITROGEN  BEHAVIOR  IN  LAND 
TREATMENT  OF  WASTEWATER 

H.M.  Selim  and  I.K.  Iskandar 


INTRODUCTION 

In  land  treatment,  nitrogen  is  almost  always  the  fac¬ 
tor  limiting  the  rate  of  wastewater  application.  Ex¬ 
cessive  nitrate  nitrogen  concentration  in  groundwater 
is  of  great  health  concern  due  to  its  association  with 
infant  methemoglobinemia  (blue  baby  syndrome)  and 
eutrophication  of  natural  waters.  Consideration  of 
land  treatment  as  an  alternative  to  advanced  waste 
treatment  has  been  hampered  by  the  lack  of  scientific 
data  on  the  fate  of  N  that  would  allow  efficient  and 
cost-effective  design  of  systems  without  incurring 
health  risks. 

The  N  behavior  in  land  treatment  is  affected  by 
numerous  physical,  chemical  and  biological  processes 
and  environmental  conditions.  Iskandar  and  Selim  (1978) 
evaluated  existing  models  for  prediction  of  NOj-N  in 
percolate  water  in  land  treatment.  They  concluded 
that  several  models  developed  to  describe  one  or  more 
processes  in  agricultural  regimes  can  be  adapted  for 
land  treatment.  However,  existing  large  models  being 
used  for  prediction  of  N  transformation  and  transport 
in  agricultural  land  must  be  modified  and  simplified  for 
use  under  land  treatment  conditions.  The  fact  that  nitro¬ 
gen  is  applied  in  small  amounts  repeatedly  (most  often 
weekly)  in  land  treatment,  in  contrast  to  normal  agri¬ 
cultural  fertilizing  practice,  should  produce  significant 
differences  in  the  nitrogen  transformation  processes. 

Also,  the  soils  under  land  treatment  are  most  often 
near  or  above  field  capacity  so  that  the  water  flow  pat¬ 
tern  as  well  as  N  transformation  processes  will  vary  sig¬ 
nificantly  from  those  of  an  agricultural  regime. 

The  objectives  of  this  report  are  to  describe  a  sim¬ 
plified  model  for  prediction  of  nitrogen  behavior  in 
slow  and  rapid  infiltration  land  treatment  systems. 


Model  development,  computer  program  listing  and 
documentation  and  sensitivity  analysis  of  model  para¬ 
meters  are  included.  Validation  of  the  developed  model 
will  be  the  subject  of  a  later  report. 

THE  MODEL 
Modeling  objectives 

The  objectives  of  developing  a  dynamic  nitrogen 
model  were  to 

1.  Develop  a  computer  model  for  N  behavior  to 
simulate  the  physical,  chemical  and  biological  processes 
in  slow  and  rapid  infiltration  systems. 

2.  Enable  prediction  of  NO3-N  concentration  in  soil 
solution  and  leachate  with  time  and  space. 

3.  Assist  in  estimating  the  application  rate  and 
schedule  of  water  and  nitrogen  to  a  land  treatment  sys¬ 
tem. 

4.  Improve  land  treatment  management  techniques 
for  better  renovation  of  waste-vater  and  less  detrimen¬ 
tal  impact  on  the  environment. 

5.  Point  out  the  area(s)  of  research  need,  based  on 
model  sensitivity  analysis  and  availability  of  information 
in  the  literature. 

Main  features 

The  main  feature  of  the  computer  program  is  that 
it  is  valid  for  uniform  as  well  as  multilayered  or  stratified 
soil  profiles.  In  addition,  the  program  is  flexible  and  is 
designed  to  incorporate  the  following  (input)  conditions 
as  desired: 

1 .  Rate  of  wastewater  application. 

2.  Duration  of  wastewater  application. 

3.  Depth  of  individual  soil  layers. 


4.  Concentration  of  ammonium  and  nitrate  in  the 
wastewater. 

5.  Wastewater  application  cycle,  i.e.  scheduling. 

6.  Soil  water  properties  and  nitrogen  transformation 
mechanisms  for  individual  soil  layers. 

7.  Plant  root  distribution  and  growth  in  the  soil. 

8.  Rate  of  nitrogen  uptake  by  plants. 

9.  Evapotranspiration  rate. 

10.  Initial  distribution  of  water  and  nitrogen  species 
in  the  soil  profile. 

General  description 

Figure  1  shows  a  block  diagram  of  the  simplified 
model  presented  in  this  report.  The  model  is  formed 
of  two  main  submodels.  The  first  is  a  water  flow  sub¬ 
model  which  describes  wastewater  infiltration,  water 
movement  in  the  soil  profile,  and  rate  of  plant  uptake 
of  water  with  soil  depth  and  time.  The  second  is  a 
nitrogen  submodel  which  describes  the  transport  and 
transformations  of  N  species  in  die  soil  as  well  as  nitro¬ 
gen  uptake  by  plants.  The  model  also  includes  several 
subroutines  which  account  for  initial  and  boundary 
conditions,  plant  root  distribution  in  the  soil,  soil  water 


properties,  and  nitrogen  transformation  processes  (ion 
exchange,  nitrification  and  denitrification).  A  detailed 
flow  chart  of  the  model  and  description  of  all  subrou¬ 
tines  are  presented  in  a  later  section. 

Water  flow  equations  and  boundary  conditions 

In  order  to  describe  the  nitrogen  transformation  and 
transport  in  saturated-unsaturated  soil  profiles  under 
transient  flow  conditions,  the  following  water  flow 
equation  (Childs  1969)  must  be  solved: 

90 /dr  =  (9/9 z)  [/C(/>)  bh/bz]-dK(h)/dz-A(z,  0) 

(D 

where 

8  =  soil  water  content  (cm3 /cm3) 
h  =  soil  water  pressure  head  (cm) 

K(h)  =  soil  hydraulic  conductivity  (cm/h'i 
A(z,  8)  =  rate  of  water  extraction  (cmJ/h  cmJ) 
t  =  time  (h) 

z  =  depth  in  the  soil  (cm). 


Figure  1  ■  Diagram  of  the  simplified  nitrogen  mode!  showing  the  water  and  the  nitrogen  submodels. 


2 


Equation  1  is  commonly  known  as  the  /j-form  of  the 
water  flow  equation.  This  equation  was  chosen  over 
the  diffusivity  form  (Selim  1978),  since  it  allows  not 
only  for  saturated-unsaturated  flow  but  it  also  allows 
soil  stratification  or  layering  of  the  soil  profile.  In 
solving  eq  1 ,  the  left-hand  term  must  be  transformed 
such  that 

bd/bt  =  (bQlbh)  (bh/bt)  =  Op  (h)  bh/bt  (2) 

where  Cap  (h)  is  the  soil  water  capacity  term  (cm'1 ) 
which  is  determined  using  the  appropriate  soil  water 
characteristic  relationship  (6  vs  h). 

In  solving  eq  1  for  multilayered  soil  profiles,  the  soil 
water  hydraulic  conductivity  K(h)  and  the  soil  water 
capacity  Cap(/i)  must  be  provided  for  each  soil  layer. 

In  addition,  the  soil  water  initial  and  boundary  con¬ 
ditions  must  be  specified.  The  initial  condition  is 
dictated  by  the  initial  distribution  of  0  or  h  in  the  soil 
profile  at  some  assumed  (starting)  time.  The  boundary 
conditions  at  the  soil  surface  and  at  some  depth  L  below 
the  soil  surface  must  be  provided. 

Soil  surface  boundary  conditions 

Two  soil  surface  boundary  conditions  are  normally 
encountered  under  field  conditions:  1 )  the  water  head 
boundary  condition  and  2)  the  water  flux  boundary  con¬ 
dition. 

Water  head  boundary  condition 

This  condition  is  used  when  water  ponding,  of  some 
height/r  above  the  soil  surface,  is  encountered.  The 
height  h  imv  be  considered  as  a  variable  with  time,  i.c. 
h(t),  in  order  to  allow  for  fluctuations  during  waste- 
water  application  and  rainfall: 

h  =  hQ(t),  at  z  =  0  (3) 

This  boundary  condition  is  also  used  when  the  soil 
surface  is  under  suction,  i.e.  the  water  content  at  the 
surface  is  below  saturation.  In  such  case,  h(t)  is  nega¬ 
tive  and  is  a  measure  of  the  soil  water  suction  (nega¬ 
tive  pressure)  at  the  soil  surface 

Water  flux  boundary  condition 

This  condition  is  imposed  when  a  constant  or  time- 
dependent  flux  (or  intensity)  q(t)  of  wastewater  (or 
rainfall)  is  applied  at  the  soil  surface.  It  also  allows  for 
evaporation  between  rainfall  or  irrigation  events.  This 
condition  can  be  written  as 


Bottom  boundary  conditions 

At  some  depth  L  below  the  soil  surface,  three  boun¬ 
dary  conditions  may  be  encountered:  1)  an  impervious 
barrier,  2)  a  soil  profile  extending  to  great  depth,  and 
3)  a  groundwater  table. 

Impervious  barrier 

rhis  boundary  condition  is  used  when  an  impermeable 
layer  (e.g.  a  heavy  clay  layer)  is  encountered  at  some 
soil  depth.  Water  flow  across  such  a  barrier  is  negligible. 
The  boundary  condition  for  an  impervious  barrier  may 
be  expressed  as 

-K (h)  bh/bz+K(h)  =  0,  at  z  =  L.  (5) 

Soil  profile  extends  to  a  great  depth 

In  this  case,  the  soil  profile  is  regarded  as  a  semi¬ 
infinite  medium.  Thus,  it  is  assumed  that  at  great  depth 
the  change  in  soil  water  suction  is  zero: 

bh/bz  =  0,  *"*•“>.  (6) 

Such  a  boundary  condition  may  be  used  if  the  soil  pro¬ 
file  is  well-drained  and  of  significant  depth. 

Groundwater  table 

If  a  water  table  is  encountered  at  some  depth  L  in  the 
soil  profile,  the  water  content  d  is  maintained  at  satura¬ 
tion  0S  at  all  times.  Therefore, 

0=0S  z  =  L,  f>0.  (7) 

Furthermore,  the  soil  water  suction  or  pressure  head  is 

h  =  0  z  -  L,  t  >  0  (8) 

In  addition  to  the  above-mentioned  boundary  con¬ 
ditions,  other  conditions  are  needed  in  order  to  de¬ 
scribe  the  water  flow  at  the  interface  between  soil  layers 
in  multilayered  or  stratified  soil  profiles.  For  example, 
we  may  consider  a  soil  profile  consisting  of  three  soil 
layers:  I,  II.  and  III.  The  length  of  each  soil  taver  is 
indicated  by  L\,  Llt  and  L3  (see  Fig.  2).  The  appro¬ 
priate  boundary  conditions  at  the  interface  between 
two  soil  layers  are  (Selim  1978). 


/  =  /.,,  t>  0, 

(9) 

/»„  =/>!„.  z  =  L-i+l-2,  t>  0, 

(10) 

where  h and  are  the  pressure  heads  in  layers 
I,  II,  and  III,  respectively.  These  boundary  conditions 


q[t)  =  -  K(h)  bhlbz+K(h),  at  /  =  0 


(4) 


rr 
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Figure  2.  Schematic  diagram 
of  a  multilayered  soil  profile. 


are  necessary  in  order  to  maintain  the  continuity  of  the 
pressure  head  h  at  the  boundary  interfaces. 

Nitrogen  transformations  and  transport  equations 
and  boundary  conditions 

In  the  development  of  this  simplified  model,  three 
major  approximations  have  been  made.  The  first 
simplification  is  that  the  nitrification  process  was  con¬ 
sidered  as  a  single  step,  i.e.  NH4  -*•  NOj  rather  than  a 
two-step  process  (NH4  -*  NOj  NOj).  Such  an  as¬ 
sumption  is  considered  adequate,  since  N02  in  most  soils 
under  neutral  pH  conditions  is  rapidly  oxidized  to  N03. 
The  second  major  simplification  is  that  the  organic-N 
phase  was  not  incorporated  in  the  model.  It  was  as¬ 
sumed  that  the  net  change  (over  a  short  period  of  time) 
in  organic-N  content  is  small  and  the  rate  of  nitrogen 
mineralization  as  well  as  immobilization  arc  extremely 
slow.  The  third  simplification  is  that  oxygen  diffusion 
in  the  soil  profile  was  not  incorporated.  Therefore, 
denitrification  of  nitrate  in  the  soil  was  assumed  to  be 
a  function  of  the  degree  of  soil  water  saturation  only. 

The  nitrogen  transformation  processes  considered 
were:  nitrification  of  NHj  to  NOj,  denitrification  of 
NOj,  and  ion  exchange  of-NHj  {seSTig.  3).  The  ion- 
exchange  process  was  assumed  to  be  instantaneous, 
whereas  nitrification  and  denitrification  processes 
were  of  the  first-order  kinetic  type  (Selim  et  al.  1976 
and  Selim  and  Iskandar  1978).  A  distribution  coeffi¬ 


cient  K0  (cm3/g)  was  used  to  describe  the  instantane¬ 
ous  (reversible)  ammonium  release  from  exchange  sites 
to  soil  solution.  The  first-order  kinetic  rate  coeffici¬ 
ents  associated  with  the  nitrification  and  denitrifica¬ 
tion  processes  were  k ,  and  k2  (h'1),  respectively. 

The  assumptions  that  these  nitrogen  transformation 
processes  follow  first-order  kinetic  reaction  were  based 
on  studies  by  McLaren  (1970, 1971 ),  Mehran  and  Tan- 
ji  (1974),  and  Hagin  and  Amberger  (1974). 

Soil  environmental  conditions  such  as  soil  suction, 
aeration,  temperature,  organic  matter  content,  and  pH 
have  significant  effects  on  the  various  nitrogen  trans¬ 
formation  mechanisms.  In  order  to  incorporate  these 
factors,  the  rate  coefficients  were  expressed  (Selim  et 
al.  1976)  as 


*7  u< 

(11) 

II 

**1 

(12) 

where  kj  and  are  considered  constants  for  each  in¬ 
dividual  soil  layer  and  f7  and  f2  are  empirical  functions 
which  describe  the  influence  of  the  previously  mentioned 
environmental  conditions  on  nitrification  and  denitri¬ 
fication,  respectively. 

The  transport  of  NH4-N  and  NOj-N  in  the  soil  solu¬ 
tion  occurs  as  a  result  of  molecular  diffusion,  mechanical 
dispersion,  and  convection  or  mass  flow.  Molecular  dif¬ 
fusion  results  from  the  random  thermal  movement  of 
molecules,  whereas  mechanical  dispersion  results  from 
the  velocity  distribution  of  water  in  the  soil  pore  space. 
For  all  soil  layers,  a  single  dispersion  coefficient  D  is 
commonly  used  which  combines  mechanical  dispersion 
and  diffusion.  Therefore,  the  convective-dispersive 
equations  governing  NH4-N  and  NOj-N  transport  may 
be  expressed  (Misra  et  al.  1974,  Davidson  et  al.  1977, 
and  Selim  and  Iskandar  1978)  as 

9(0C)/9r  =  (9/9z)  (6D  9C/9z)-9(vC)/9z 


-6  k  j  C-p  dS/dt-q  nh4  (13) 

d(6Y)lbt  =  (9 /dz)  ( OD  9K/9z)-9(vV')/9z 

H)k  jC-dk2Y-q^Q3  (14) 


where  C  =  concentration  of  NH4-N  in 
soil  solution  Oug/cm  ) 

Y  -  concentration  of  NO3-N  in 
soil  solution  (pg/cm3) 
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Figure  3.  Schematic  diagram  of  the  nitrogen  transformation  processes 
considered  in  the  nitrogen  submodel. 


D  =  solute  dispersion  coefficient 
(cm2/h) 

v  =  soil  water  flux  (cm/h) 

S  =  amount  of  NH4  in  the  exchange¬ 
able  phase  per  gram  of  soil 

(fM 

p  =  soil  bulk  density  (g/cm3) 
kj  and  k2  -  kinetic  rate  coefficients  for 

nitrification  and  denitrification 
(  'M.  respectively 

<7NH4  and  pN03  =  raus  of  plant  uptake  of  NH4-N 
and  NO3-N  per  unit  soil  vol¬ 
ume  (pg/cm3  h),  respectively. 

The  first  two  terms  on  the  right-hand  side  of  eq  13 
and  14  account  for  solute  transport,  and  are  usually 
called  the  dispersion  and  mass  flow  terms,  respectively. 
The  third  and  fourth  terms  of  eq  13  account  for  nitri¬ 
fication  and  ion  exchange,  respectively,  of  NH4-N. 
Similarly,  the  third  and  fourth  terms  of  eq  14  represent 
the  nitrification  and  denitrification  processes,  respec¬ 
tively.  The  ion  exchange  process  governing  NH4-N 
adsorption-desorption  was  assumed  to  be  of  the  linear 
Freundlich  type,  i.e. 

S  =  KDC,  or  3 S/dt  =  KD  bC/bt,  (1 5) 

where  KD,  commonly  called  the  distribution  coefficient 
(cm3/g),  represents  the  ratio  between  the  amount  of 
NH4-N  adsorbed  and  its  concentration  in  the  soil  solu¬ 
tion. 

Rearrangements  of  ea  13  and  incorooration  of  eq 
1 3  yield  the  following  simplified  equation  for  am¬ 
monium  transport  and  transformation: 


R  bC/bt  =  D  b2C/bz2-(V/6)  bC/bz-k,  C-(qNm/0) 

(16) 

where  R  is  the  retardation  factor  for  ammonium  ex¬ 
change: 

R=  ]  +  pKD/e.  (17a) 

and  V  is  expressed  as 

V  -  v-D  bO/bz.  (17b) 

Similarly,  eq  14  after  rearrangements  yields  the  following 
equation  for  nitrate  transport  and  transformations: 

bYlbt  =  D  b2Ylbz2-(V/e )  bYlbz+k,C-k2Y-(qNO3l0). 

(18) 

It  should  be  noted  that  in  the  case  of  multilayered 
soil  profiles,  soil  water  and  nitrogen  transformation 
parameters  (e.g.  p,  K(h),  kp,  kltk2,  etc.)  must  be  pro¬ 
vided  for  each  individual  soil  layer. 

To  solve  the  ammonium  and  nitrate  transport  and 
transformation  equations,  eq  16  and  18,  the  initial  and 
boundary  conditions  must  be  specified.  During  waste- 
water  application,  the  soil  surface  boundary  conditions 
for  eq  16  and  18,  respectively,  are: 

vCs  =  -4DbClbz+vC,  z  =  0,  t  <  T,  (19) 

and 


vYs  =  -6D  b  Y/bz+v Y,  z  =  0,  t  <  T,  (20) 

where 

Cs  and  Ys  =  NH4-N  and  NO3-N  concentrations 
in  applied  wastewater,  respectively 
(fig  N/ml) 

v  =  q(t),  flux  or  intensity  of  wastewater 
application  (cm/h) 

T  =  duration  of  wastewater  application 

(h). 

The  above  equations,  eq  19  and  20,  are  commonly 
called  the  dispersive-convective  boundary  conditions 
which  can  also  be  applied  to  describe  rainfall  events. 

In  such  a  case,  the  intensity  and  duration  of  rainfall 
(v  and  T)  must  be  specified  and  Cs  and  Ys  in  rainwater 
may  be  considered  zero. 

Following  the  termination  of  a  wastewater  applica¬ 
tion  or  rainfall  event  (t  >  T),  the  surface  boundary  con¬ 
ditions  become 

bC/bz  =  0,  z=0tt>T  (21 ) 

and 

bY/bz  =  0,  z  =  0,t>T  (22) 

Furthermore,  the  boundary  condition  at  the  bottom 
of  the  soil  profile  (z  =  L)  is 

bC/bz  =  0,  z  =  L,t>  0  (23) 

bY/bz  =  0,  z  =  L,t  >  0  (24) 

In  addition,  the  boundary  conditions  at  the  boundary 
interface  between  two  soil  layers  (see  Fig.  2)  may  be 
written  as 

C|  =  C|t  and  K|  =  y(|,  z-Ly,t>  0  (25) 

^ll  =  ^lil  anc^  ^ll  =  ^lll'  z-Li+L.2,t>0 

(26) 

where  the  subscripts  I,  II,  and  III  refer  to  the  first,  se¬ 
cond.  and  third  soil  layers,  respectively.  Similar  to 
the  equations  for  water  flow,  eq  25  and  25  above 
are  needed  in  order  to  maintain  the  continuity  of 
NFI4-N  and  NO3-N  concentrations  at  the  boundary 
interface. 


Water  and  nitrogen  uptake  by  plants 

Plant  uptake  of  water  and  nitrogen  from  the  soil 
root  zone  is  an  important  factor  in  the  renovation  of 
wastewater  applied  to  soil.  Recent  studies  have  shown 
that  in  a  slow  infiltration  land  treatment  system,  a 
major  portion  of  applied  wastewater  nitrogen  (up  to 
70%)  was  taken  up  by  plants  (Iskandar  et  al.  1976). 
Therefore  in  modeling  the  fate  of  nitrogen  in  soil,  it  is 
important  to  incorporate  a  plant  uptake  model  that 
provides  accurate  predictions  of  the  rate  of  plant  uptake 
during  the  growing  season.  However,  as  Nye  and  Tinker 
(1977)  pointed  out,  the  major  difficulties  in  modeling 
of  plant  uptake  are  the  lack  of  quantitative  measurements 
on  the  root  development  and  distribution  as  well  as  the 
inaccuracy  of  soil  physical  measurements. 

At  present  there  are  two  approaches  for  modeling 
plant  root  uptake  of  water  and  nutrients  in  soils:  1 ) 
a  “microscopic”  approach  where  the  water  and  nutrient 
flux  to  a  single  root  is  considered  (Nye  and  Marriot 
1969,  Claassen  and  Barber  1976),  and  2)  a  “macro¬ 
scopic”  approach  where  the  root  system  as  a  whole  is 
considered  (Molz  and  Remson  1970,  Davidson  et  al. 
1977,  Selim  and  Iskandar  1978).  In  this  simplified 
model  the  macroscopic  approach  is  used  to  describe  the 
water  as  well  as  the  nitrogen  uptake  by  plant  roots.  The 
extraction  or  sink  term  A(z,  0)  for  water  uptake  (eq  1) 
is  represented  as 

Z 

A(z,6)  =  TR(z)K(h)l  £  R(z)K(h)dz  (27) 

where  Z  is  the  maximum  depth  of  the  root  zone  in  the 
soil  (cm)  and  T  the  evapotranspiration  rate  per  unit  area 
of  soil  surface  (cm/h).  The  term  R(z)  is  the  root  dis¬ 
tribution  as  a  function  of  depth  in  the  soil  profile.  Spe¬ 
cifically  the  root  distribution  R(z)  is  the  length  of  roots 
(cm)  as  a  function  of  soil  depth  and  time.  Equation  27 
was  proposed  by  Molz  and  Remson  (1970)  and  was 
successfully  used  in  predicting  the  water  uptake  when 
the  evapotranspiration  rate  T  was  met.  Such  conditions 
are  satisfied  when  high  soil  water  contents  (low  suctions) 
are  maintained  in  the  soil  root  zone,  such  as  in  land 
treatment-slow  infiltration  systems. 

The  terms  <?NH4  anc*  9no3  in  e9  16  and  18  acc<  unt 
for  the  rate  of  uptake  of  NH4-N  and  NOj-N,  respec¬ 
tively.  Here  the  Michaelis-Menten  approach  was  used 
to  determine  the  rate  of  N  uptake  as  a  function  of 
root  density  and  concentration  of  ammonium  and  ni¬ 
trate  in  the  soil  solution.  Therefore  the  rate  of  N  up¬ 
take  may  be  expressed  as: 
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<?NH4  =  W  Cl[Km+{C+Y)\ 

(28) 

0NO3='max  W*m+(C+V1] 

(29) 

In  eq  28  and  29,  /max  is  the  maximum  rate  of  N 
uptake  per  unit  root  length  (/ig/h  cm)  when  the  con¬ 
centration  of  nitrogen  in  the  soil  solution  is  extremely 
high,  and  the  term  Km  is  the  Michaelis  constant  (jug/ml) 
which  is  the  concentration  of  N  at  Vi  /max.  Both  /max 
and  Km  are  determined  by  measuring  N  uptake  in 
solution  cultures  having  different  nitrogen  concentra¬ 
tions  (daassen  and  Barber  1976).  In  this  model  the 
values  of  /max  and  Km  were  considered  similar  for  both 
ammonium  and  nitrate  uptake*. 


NH4-N  transport  and  transformation  (eq  16)  may  be 
expressed  as 

Rn+1  [Cn+1_cn]  =  yQ  ^1  _2C|*1  ) 

no[c^1-2c|’^1) 

-(we>r+l  flcsv-cr’i 

-Atk,C?-At(qm4/0)?  (31) 

and  the  finite  difference  approximation  for  N03-N  (eq 
18)  is 


Method  of  model  solution 

The  water  and  nitrogen  equations  (eq  1, 16  and  18) 
are  nonlinear  partial  differential  equations  and  cannot 
be  solved  analytically.  Therefore,  these  equations,  sub¬ 
ject  to  the  above  described  initial  and  boundary  condi¬ 
tions,  were  solved  using  numerical  analysis  techniques. 
The  method  of  solution  v  as  by  explicit-implicit  finite 
difference  approximation  \nenrici  1962,  Varga  1962, 
and  Carnahan  et  al.  1 969).  i  ms  meuiod  was  success¬ 
fully  used  by  Selim  (1978)  for  transient  water  and 
solute  movement  in  multilayered  soil  profiles.  Finite 
difference  approximations  provide  distributions  of 
soil  water  content,  water  suction,  NH4-N  and  NOj-N 
concentrations  at  incremental  distances  A z  in  the  soil 
profile,  and  discrete  time  steps  At.  In  finite  difference 
form,  a  variable  such  as  A  is  expressed  as  A"  =  h(z,  t ) 

=  h(iAz,  nAt)  where  /  =  1, 2,  3, ...,  /,  and  n  =  1,2,... 
Therefore,  the  finite  difference  approximation  for  the 
water  flow  equation  (eq  1 )  is 

Cap  (/»P+,/2)  [/i"+1-/7"l  =yK(h£$)  [A^r1 1 

+T*«lJ£)  K,-AP) 

-|3[^(A  J.,^)  +Ar(AP  1/2)|  -At  A  "  (30; 

where  y  -  At/2{Az)2  and  (3  =  At/Az.  Similarly,  finite 
difference  approximations  for  the  equation  governing 


*S.A.  Barber,  Department  of  Agronomy,  Purdue  University, 
personal  communication  1979. 


VP+1  -  Vi1  =  y D  |  -2  V'P+1  +VP_V  : 


no(v^1-2fP+vPt] 

-(l//0)P+1  (3[V"Y-KP+M 

+Af*,  CP-ArA?  VP-Ar(qNO3/0)P.  (32) 

Equations  30,  31 ,  and  32  are  nonlinear  since  Cap 
(A"+1 )  and  K(AP+,/2)  are  dependent  on  h"+^2  for  which 
solutions  are  being  sought.  The  iteration  method  de¬ 
scribed  bv  Remson  etal.  (1971)  is  usually  used  to  pre¬ 
dict  An+1'2  usingAn.  Selim  and  Kirkham  (1973)  showed 
that  solution  of  the  water  flow  equation  can  be  approxi¬ 
mated  satisfactorily  using  An  when  smaller  values  of 
At  than  required  for  stable  solution  are  used.  This 
simplifies  the  computation  considerably  since  the  sys¬ 
tem  of  equations  becomes  linear.  Accordingly,  the 
approximations,  K(h"+]/2)  =  AC(AP)  and  Cap  (AP+^2) 

=  Cap  (A"),  were  made. 

Incorporation  of  initial  and  boundary  conditions  ir 
their  finite  difference  forms  and  rearrangement  of  eq 
30,  31 ,  and  32  yield  three  linear  systems  of  equations. 
Rearranging  the  finite  difference  form  of  water  flow 
equation  (eq  30)  yields 

«VArV  Cl1  (33) 

where  rfP  =  -y  AC(AP)/2) 

eP  =  Cap  (AP)n  [K(Aki/2)+Af(A?_,/2).l 

^P  =  -7  K(h^.y2) 

wP  =  Cap(AP)A"nK(AP1/2)AP, 

-7  [K(A^/2)+/C<AnM/2)]  AJ1 
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-PlKKv2)+K(hly2)]-*tA?. 


By  including  the  initial  and  boundary  conditions  in 
their  finite  difference  forms,  eq  33  can  be  written  in 
matrix-vector  notation  as 


Bh^' =w  (34) 

where  B  is  a  tridiagonal  real  matrix  and  h  and  w  denote 
the  associated  real  column  vectors  (the  arrows  indicate 
vectors).  The  matrix  B  may  be  written  as 


„n  _n 
*1*1 

<*2*2*2 

<*3*3*3 


B  = 


„n  „n 
*1-1  *  1-1 


d 


n 

I 


(35) 


The  coefficients  of  the  main  diagonal  of  the  matrix  B, 
in  absolute  values,  are  greater  than  the  raw  sum  of  the 
off-diagonal  coefficients.  Hence,  the  matrix  B  is  strict¬ 
ly  diagonally  dominant  (Varga  1962,  p.  23).  Therefore 
the  matrix  is  nonsingular,  and  there  exists  a  solution 
h  n+1  for  the  matrix  vector  equation  (eq  34)  that  is 
unique.  The  tridiagonal  system  of  equations  was  solved 
by  an  adaptation  of  the  Gaussian  algorithm  as  described 
by  Henrici  (1962,  p.  352).  The  second  and  third  sys¬ 
tems  of  eq  31  and  32  for  NH4-N  and  NOj-N  transport 
and  transformation,  respectively,  were  solved  similarly. 

It  is  obvious  from  eq  16  and  18  that  nitrogen  trans¬ 
port  is  dependent  on  0  and  v,  both  of  which  are  variable 
under  transient  flow  conditions.  Thus,  in  order  to  de¬ 
scribe  nitrogen  transport  for  transient  conditions,  the 
water  flow  equation  (eq  1 )  must  be  solved  prior  to  the 
nitrogen  transformation  and  transport  equations.  There 
fore,  for  any  time  step,  n+1 ,  where  all  variables  are 
assumed  to  be  known  at  time  step  n,  eq  30,  31 ,  and  32 
are  solved  sequentially  until  a  desired  time  t  is  reached. 


MODEL  SENSITIVITY 

In  order  to  provide  a  complete  sensitivity  analysis 
of  model  parameters,  it  is  essential  to  investigate  each 
parameter  separately.  This  is  usually  achieved  by  first 
studying  the  influence  of  each  parameter  on  model 
results  for  a  wide  range  of  values,  with  all  other  para¬ 
meters  remaining  unchanged.  Second,  when  two  or 


more  parameters  prove  to  be  more  significant  compared 
to  other  model  parameters,  such  two  or  more  parameters 
are  investigated  simultaneously  for  a  range  of  values. 

For  the  model  presented  here  a  complete  sensitivity 
analysis  was  not  attempted  since  a  large  number  of 
model  parameters  were  involved.  Therefore,  it  was 
decided  to  study  only  selected  parameters  which  were 
chosen  for  sensitivity  analysis.  These  parameters  are 
1 )  the  rate  of  nitrification  (ky ),  2)  the  distribution 
coefficient  (KD)  for  NH4-N  ion-exchange,  and  3)  rate 
of  plant  uptake  (/max).  The  influence  of  different 
nitrogen  transformation  parameters  for  individual  soil 
layers  and  the  rate  of  denitrification  were  not  investi¬ 
gated.  Moreover,  soil  water  properties,  schedule  of 
wastewater  application  and  NH4-N  and  NO3-N  con¬ 
centrations  in  the  wastewater  remained  unchanged  for 
all  cases  studied. 

The  simultaneous  transport,  transformation,  and 
plant  uptake  of  nitrogen  and  water  were  simulated  using 
the  model  presented  here.  Input  soil  and  water  para¬ 
meters,  initial  and  boundary  conditions,  etc.,  were 
similar  to  those  of  the  CRREL  research  experimental 
facilities.  Wastewater  was  assumed  to  contain  25  pg/ml 
of  NH4-N  and  zero  nitrate  content.  The  schedule  of 
application  was  0.5  or  1  week  and  the  total  amount  of 
wastewater  applied  was  5  cm  per  each  application.  The 
soil  parameters  chosen  were  for  a  Windsor  sandy  loam 
soil  having  three  distinct  soil  layers.  Chemical  charac¬ 
teristics  of  this  soil  are  presented  elsewhere  (Iskandar 
et  al.  1979,  Iskandar  and  Nakano  1978).  The  total 
length  of  the  soil  profile  was  assumed  to  be  1  50  cm  and 
the  thicknesses  of  the  first,  second  and  third  layers  were 
1 5,  30  and  1 05  cm,  respectively.  A  water  table  (zero 
water  pressure  head)  was  assumed  at  the  bottom  (1 50 
cm  depth)  of  the  soil  profile.  The  soil  water  properties 
for  each  soil  layer  were  described  by  the  following 
mathematical  expressions: 


0(h)  =  OJ\U(-hla)b\ 

(36) 

K(0)  =  rj  exp  (a0). 

(37) 

The  parameters  for  each  soil  layer  were  obtained  by 
"best  fit”  of  experimental  data  (Iskandar  and  Nakano 
1978).  The  values  of  a  and  b  were  100  and  1  for  the 
first  layer,  40  and  1  for  the  second  layer,  and  30  and 
1  for  the  third  layer,  respectively.  The  values  for  para¬ 
meters  r?  and  a  were  9.6x  10‘4  and  27.63  for  the  first 
layer,  2.2x  10'4  and  30.7  for  the  second  layer,  and  2.1 
x104  and  38.87  for  the  third  layer,  respectively.  The 
values  for  water  content  at  saturation  0S  were  0.44, 
0.42  and  0.34  (cm3/cm3)  for  the  first,  second  and 
third  layers,  respectively.  Furthermore,  the  values  for 
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soil  bulk  density  p  were  1 .41  for  the  first  layer,  1 .59 
for  the  second  layer,  and  1.55  g/cm3  for  the  third  layer. 

Several  nitrogen  transformation  rate  coefficients 
and  uptake  rates  were  chosen  in  order  to  illustrate  the 
significance  of  these  processes  on  the  fate  of  wastewater 
nitrogen  in  land  treatment  systems.  The  rate  of  nitri¬ 
fication  ranged  from  0.025  to  0.5  h'1  whereas /t2 
was  maintained  constant  at  0.01  h’1 .  For  the  nitrifica¬ 
tion  kinetic  reaction,  the  function  f,  (eq  1 )  which  de¬ 
scribes  the  dependence  of  the  reaction  on  soil  environ¬ 
mental  conditions  remained  unchanged  and  was  ex¬ 
pressed  as  a  function  of  pressure  head  (Hagin  and  Am- 
berger  1974): 


'l  = 


0 

,0.005  (-*- 10) 
'0.2+0.006  (-Tr-40) 
l0.5+0.0015  (-*-100) 
1.0-0.002  (-*-433) 


*  >  -10  cm 

*  >  -50  cm 
h  >  -1 00  cm 
h  >  -433  cm 
h  <  -433  cm. 


The  denitrification  function  f2  (eq  12)  was  considered 
as  a  function  of  the  degree  of  water  saturation  in  the 
soil: 


f2  =  0  for  (0/0s)  <  0.8 

f2  =  (0-0.8  0S)/O.1  0S  for  0.8  <  0/05  <  0.9 
f2  =  1  for  0/0,  >  0.9 


Soil-Water  Sue  non  (cm) 


Figure  4.  Simulated  soil  water  suction  distributions 
in  a  Windsor  soil  for  one  week  following  applica¬ 
tion  of  5  cm  of  wastewater 


where  0S  is  soil  water  content  at  saturation  (cm3/cm3). 
Furthermore,  the  value  of  the  coefficient  KD  for  the 
ammonium  exchange  ranged  from  0.05  to  1 .5  cm3/g. 
All  these  parameters,  for  the  cases  shown  here,  were 
similar  for  all  three  soil  layers. 

The  root  distribution  R{z)  used  in  the  sensitivity 
analysis  was  in  the  form  of  an  exponential  expression: 

R{z)  =  226  exp  (-0.1  z) 

This  exponential  function  provides  a  sharp  decrease  of 
the  root  distribution  (or  density)  with  soil  depth  z.  In 
this  case,  77.8%  of  the  roots  were  actually  observed  to 
be  in  the  top  1 5  cm,  1 1.2%  in  the  1 5-  to  30-cm  soil 
depth  and  only  5%  of  the  roots  in  the  30-  to  60-cm  soil 
depth.  This  mathematical  expression  for  root  distribu¬ 
tion  was  obtained  from  field  measurements  of  root 
length  with  depth  (Iskandar  unpublished  data  1979). 

It  represents  a  two-year-old  mixture  of  reed  canary- 
grass  and  tali  fescue  grass  irrigated  with  5  cm  of  waste- 
water  weekly.  Furthermore,  a  value  of  1 .0  ppm  for 
the  Michaelis  constant  Km  was  chosen  for  all  cases 
presented  here.  Values  of  lmix  ranged  from  0.0005 
to  0.001 5  ftg  N/h  cm  of  roots. 

In  order  to  illustrate  the  capabilities  of  this  model, 
simulated  water  content  and  nitrogen  distributions  in 


the  soil  are  presented  in  Figures  4,  5  and  6  for  one  week 
following  the  application  of  5  cm  of  wastewater.  In 
this  case  the  rate  coefficients  for  nitrification  were 
taken  as  0.1  h'1  and  KD  as  0.25  cm3/g.  Following  the 
application  of  5  cm  of  wastewater,  which  contained 
25  ppm  of  NH4-N,  we  see  that  NH4-N  concentration 
in  the  soil  solution  decreased  drastically  with  time.  This 
disappearance  of  NH4-N  was  due  to  its  conversion  to 
NO3-N  through  nitrification  and  to  uptake  by  plants. 
The  NO3-N  distributions  shown  in  Figure  6  show  that 
the  maximum  peak  concentration  occurred  on  day  1 . 
This  is  primarily  due  to  the  initial  NOj-N  in  the  soil 
profile.  Furthermore,  the  initial  peak,  which  was  located 
at  the  1 5-cm  depth,  penetrated  to  a  depth  of  35  cm 
during  the  first  day.  Downward  movements  of  the 
nitrate  peak  to  lower  soil  depths  continued  with  time 
at  a  decreasing  rate.  This  slow  movement  of  nitrate 
in  the  soil  profile  was  due  to  the  continuing  decrease 
of  soil  water  flow  during  water  redistribution,  i.e. 
following  the  infiltration  of  applied  wastewater.  Mean¬ 
while,  the  nitrate  concentration  of  the  leachate  (z  =  1 50 
cm)  continued  to  increase  with  time. 

Figure  7  shows  the  cumulative  nitrogen  uptake  by 
plants  with  time  during  weekly  application  of  waste- 
water.  Simulated  results  show  that  NH4-N  uptake  was 
much  greater  than  that  for  N03-N  at  all  times.  This 
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I  max  =  0.001  pg  N/h  cm.  Arrows  indicate  wastewater  appli¬ 
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Figure  8.  Simulated  concentration  distributions  ot 
N03-N  in  the  soil  profile  where  weekly  application 
of  5  cm  of  wastewater  was  maintained.  Model 
parameters  were  k,  =  0.025  h~l,  K n  =  0.25  cm3  /g 
and  I  max  =  0.001  pg  N/h  cm. 


Figure  9.  Simulated  concentration  distributions  ot 
NOi-N  in  the  soil  profile  where  weekly  application 
of  5  cm  of  wastewater  was  maintained.  Model 
parameters  were  k,  =  0.05  h'1,  Kp  =  0.25  cm3/g 
and  I  max  =  0. 001  pg  N/h  cm. 


is  due  to  the  low  concentration  of  N03-N  in  the  top 
portion  of  the  soil  profile  for  an  extended  period  of 
time.  In  contrast,  the  NH4-N  remained  at  shallow  soil 
depth  due  to  ion-exchange  which  resulted  in  higher 
NH4-N  uptake.  It  should  be  noted  that  in  the  simu¬ 
lation,  77.8%  of  plant  roots  were  in  the  top  1 5  cm  of 
the  soil  profile. 

Figures  8-1 2  show  the  influence  of  changing  the 
nitrification  rate  coefficients*,  on  the  nitrate  distribu¬ 
tion  in  the  soil  profile  during  weekly  application  of 
wastewater.  In  the  cases  considered,  the  values  of  k , 
range  from  0.025  to  0.5  h'1 ,  whereas  KD  and  lmix  re¬ 
mained  constant  (KD  =  0.25  cm3/g  and  lmax  =  0.001 
pg-N/h  cm).  The  simulated  results  show  that  as  the 
rate  of  nitrification  increased,  the  maximum  peak  con¬ 
centration  also  increased  at  the  end  of  the  weekly 
wastewater  applications.  These  increases  in  peak  con¬ 
centrations  were  obviously  due  to  the  imposed  increase 
in  faster  conversion  of  NH4-N  to  N03-N  in  the  soil. 
Furthermore,  maximum  peak  concentrations  were  locatec 
at  shallower  soil  depths  as*,  increased.  The  location 
of  the  peak  closer  to  the  soil  surface  coincides  with  that 
for  NH4-N.  The  simulated  results  also  show  that  the 
rate  of  nitrification  directly  influences  the  effluent  con¬ 
centration  (at  /  =  1 50  cm).  As  expected,  the  effluent 
NOj-N  concentration  (at  any  particular  time)  increased 


ar  the  rate  of  nitrification  increased. 

Figure  1 3  shows  the  effect  of  changing  the  rate  of 
nitrification  on  the  total  amounts  of  NH4-N  and 
NOj-N  in  the  soil  profile  during  weekly  wastewater  ap¬ 
plication.  These  simulated  results  were  obtained  by  in¬ 
tegrating  the  nitrogen  concentration  profiles  shown  in 
Figures  8-1 2.  As  expected,  as  the  rate  of  NH4-N  con¬ 
version  to  NOj-N  increased,  the  total  amount  of  NH4-N 
in  the  soil  profile  decreased  and  the  total  amount  of 
NOj-N  in  the  soil  profile  increased.  It  should  also  be 
emphasized  here  that  the  relationship  between  N  in 
the  soil  profile  and  *,  (Fig.  13)  is  a  nonlinear  one.  This 
is  clearly  shown  by  the  leveling  of  the  curves  for  *, 
greater  than  0.2  h'1 .  Such  leveling  of  the  curves  indi¬ 
cates  that,  for  the  cases  considered,  the  influence  of 
*,  on  the  fate  of  nitrogen  is  negligible  for  *,  >  0.2  h*1 . 
Therefore,  it  is  expected  that  beyond  such  *,  values, 
other  factors  which  affect  the  amount  of  NH4-N  in  the 
soil  profile  (such  as  the  rate  of  N  uptake,  the  ion-ex¬ 
change  of  NH4-N,  and  the  amount  of  weekly  waste- 
water  applied)  will  influence  the  NO3-N  in  the  soil  pro¬ 
file.  For  example,  if  the  amount  of  weekly  wastewater 
were  50  rather  than  25  ppm,  it  would  be  reasonable  to 
expect  that  the  leveling  will  be  at  a  higher  *,  value 
than  0.2  h"' .  The  opposite  would  be  true  for  a  reduced 
rate  of  plant  uptake. 
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Figure  14.  Simulated  concentration  distributions  of 
NOyN  in  the  soil  profile  where  weekly  application 
of  5  cm  of  wastewater  was  maintained.  Model 
parameters  were  k,  =  0.1  h~' ,  KD  =  0.05  cmi/g 
and  lmox  =  0.00 1  pg  N/h  cm. 

The  influence  ot  a  range  or  K0  values  on  the  fate 
of  wastewater  nitrogen  is  shown  in  Figures  14-16.  The 
KD  values  considered  ranged  from  0.05  to  1.5  cm3/g 
which  simulates  a  wide  range  of  soils  having  different 
cation  exchange  capacities.  The  kt  and  /max  values 
remained  constant  (fc,  =  0.1  h'1  and  =  0.001 
pg-N/cm  h)  for  these  cases.  For  small  values  of  KD 
(Fig.  14),  the  nitrate  concentration  profiles  showed  a 
rapid  movement  in  the  soil  profile.  Such  rapid  nitrate 
leaching  is  a  direct  result  of  low  NH4-N  retardation  in 
the  soil  profile.  In  contrast,  for  larger  KD  values  (Fig. 

1 5  and  1 6),  the  retardation  of  NH4-N  to  transport  in 
the  soil  profile  was  greater,  resulting  in  slower  move¬ 
ment  of  NOj-N  in  the  soil  profile.  This  slow  move¬ 
ment  of  nitrate  nitrogen  is  illustrated  by  the  location 
of  the  peaks  at  shallow  soil  depths  as  well  as  higher 
peak  concentrations. 

The  effect  of  lmax,  the  maximum  rate  of  N  uptake, 
on  the  fate  of  nitrogen  in  the  soil  profile  as  well  as  the 
cumulative  N  uptake  are  shown  in  Figures  1 7-19.  Here 
the  range  of  ^  was  from  0.005  to  0.001 5  pg-N/n  cm 
of  root  length.  All  other  parameters  were  constant: 

=0.1  h'(  and  KD  =0.25cm3/g.  As  expected,  for 
small  values  of  (Fig.  1 7),  the  N03-N  distributions 
show  significantly  higher  concentrations  than  for  large 
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Figure  15.  Simulated  concentration  distributions  o> 
NOi-N  in  the  soil  profile  where  weekly  application 
of  5  cm  of  wastewater  was  maintained.  Model 
parameters  were  k,  =0.1  h~' ,  K0  =  0.5  cm3/g 
and  lmox  =  0.001  pg  N/h  cm. 

^max  values.  Differences  in  concentration  distributions 
are  clearly  shown  in  the  upper  portion  of  the  soil  pro¬ 
file.  For  I  max  =  0.001 5  pg  N/h  cm,  plant  uptake  resulted 
in  considerable  depletion  of  the  N03-N  in  the  root  zone 
(0-30  cm).  In  contrast,  considerable  amounts  of  N03-N 
remained  in  the  top  soil  profile  for  the  case  /max  = 
0.0005  pg-N/h  cm.  The  cumulative  nitrogen  uptake  pat¬ 
terns  for  different  values  (Fig.  19)  show  that  for  /max  = 
0.005  the  uptake  with  time  was  linear.  Such  a  linear 
relationship  indicates  a  constant  rate  of  nitrogen  up¬ 
take  with  time.  For  this  case,  the  maximum  rate  of 
nitrogen  uptake  was  met,  indicating  an  abundance  of 
nitrogen  in  the  soil  root  zone  at  all  times.  However, 
for  large  values  of  the  amount  of  nitrogen  in  the 
root  zone  was  limited  as  indicated  by  the  nonlinear  up¬ 
take  patterns  for  of  0.001  and  0.001 5  pg-N/h  cm. 
Immediately  following  each  weekly  application,  the 
simulated  results  also  show  a  rapid  increase  in  the  nitro¬ 
gen  uptake  as  a  result  of  the  newly  added  nitrogen  in 
the  wastewater.  However,  three  to  four  days  after  ap¬ 
plication  of  wastewater,  when  the  nitrogen  concentra¬ 
tion  in  the  root  zone  becomes  small,  the  rate  of  up¬ 
take  decreases  drastically  with  time.  This  change  of  the 
rate  of  uptake  with  time  is  clearly  illustrated  in  Figure 
19  by  the  sudden  increase  followed  by  a  leveling  of  the 
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Figure  1 6.  Simulated  concentration  distributions  of 
N03-N  in  the  soil  profile  where  weekly  application 
of  5  cm  of  wastewater  was  maintained.  Model 
parameters  were  k,  =0.1  h  ' ,  Kc  =  7.5  cm3/g 
and  \max  =  0.001  pg  N/h  cm. 
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Figure  7  7.  Simulated  concentration  distribution  of 
NurN  in  tne  soil  profile  where  k,  =  0.1  h  ‘  ,KD  - 

0.25  cmVg  and  \mox  =  U.OUU5  pg  N/h  cm. 
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Figure  18.  Simulated  concentration  distribution  of 
NOj-N  in  the  soil  profile  where  k}  =  0.1  h'\  KD= 
0.25  cm3 /g  and  \max  =  0.0015  pg  N/h  cm. 


Figure  19.  Cumulative  N  uptake  with  time  for  \max  of 0.0005, 
0.0010 and  0.0015  pg  N/h  cm. 
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Figure  20.  Simulated  concentration  distributions  of 
NOrN  in  the  soil  profile  where  5  cm  of  wastewater 
was  applied  every  A  week.  Model  parameters  were 
the  same  as  those  used  for  Figure  10. 

cumulative  nitrogen  profile  during  each  application  ot 
wastewater.  Obviously  maximum  nitrogen  uptake  was 
not  met  for  the  cases  with  of  0.001  and  0.001 5 
ug  N/h  cm 

Figures  20  and  21  snow  the  influence  ot  the  amounts 
and  scheduling  of  wastewater  application  on  nitrate  dis 
tributions  in  the  soil  profile.  Figure  20  shows  N03-N 
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Figure  21.  Same  as  Figure  20  except  for  weekly 
application  of  10  cm  of  wastewater. 

distribution  in  the  soil  where  5  cm  ot  wastewater  was 
applied  every  3.5  days,  i.e.  twice  weekly.  This  rate  is 
twice  the  application  rate  considered  in  previous  cases 
where  only  a  5-rm  wastewater  application  was  main¬ 
tained  every  week.  Figure  21  is  for  a  different  case 
where  10  cm  of  wastewater  is  applied  in  one  application 
every  week.  All  other  parameters,  for  the  purpose  of 
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Figure  22.  Simulated  concentration  distributions  of 
N03-N  in  the  soil  prof  He  where  weekly  application 
of  5  cm  of  wastewater  was  maintained.  Model 
parameters  were  the  same  as  those  used  for  Figure 
10. 


Figure  23.  Cumulative  N  uptaRe  with  time  where  weekly 
application  of  5  cm  of  wastewater  was  maintained  for  a 
period  of  9  weeks.  Dashed  line  is  for  a  case  where  5  cm 
of  water  (or  simulated  rainfall),  rather  than  wastewater, 
was  applied  on  the  fifth  week  (day  28). 


comparison,  are  the  same  as  those  used  for  Figure  10. 
Comnarisnn  of  Figures  10  and  20  shows  that  increasing 
the  schedule  of  application  to  5  cm/'/i  week,  rather 
than  one  week,  resulted  in  more  leaching  of  the  nitrate 
nitrogen  from  the  soil  profile.  In  addition,  peak  con¬ 
centrations  were  at  lower  soil  depths  when  the  waste- 
water  was  applied  twice  weekly.  Comparison  of  Fig¬ 
ures  20  and  21  shows  that  the  total  NO3-N  in  the  soil 
was  applied  once  every  week  than  for  5  cm  every  Vi 
week.  Increased  total  NO3-N  in  the  soil  profile  (shown 
in  Fig.  21)  was  probably  due  to  increased  depth  of 
penetration  of  NO3-N  in  the  soil  for  the  10  cm/weck 
application  rate. 

rhe  influence  of  rainfall  events  on  NO3-N  distribu¬ 
tions  in  the  soil  profile  is  shown  in  Figures  22  and  23. 
Flcre,  we  compare  a  case  where  5  cm  of  wastewater  is 
applied  weekly  for  a  total  period  of  9  weeks  to  another 
case  where  in  the  fifth  week  (day  28)  5  cm  of  water  or 
simulated  rainfall  containing  no  NO3-N  or  NFI4  N  is 
applied.  Figure  22  shows  that  the  results  of  NO}-N 
distributions  in  the  soil  for  the  two  cases  were  identical 


for  the  time  of  simulation  considered.  This  was  not 
surprising  since  the  lop  portion  of  the  soil  profile  con¬ 
tained  very  little  NO3-N  at  the  end  of  the  fourth  week. 
As  a  result,  the  effect  of  the  simulated  rainfall  event  on 
the  fifth  week  was  in  limiting  the  amount  of  N  uptake 
during  that  time.  This  is  clearly  illustrated  in  Figure  23, 
where  after  the  fifth  weet  the  cumulative  N  uptake  was 
approximately  1 25  kg  N/ba  higher  for  the  case  where 
weekly  wastewater  was  maintained.  Moreover,  there 
was  a  leveling  of  N  uptake  during  the  fifth  week  for 
the  case  where  5  cm  ot  water  (or  rainfall)  was  applied 
on  day  28.  This  difference  in  the  cumulative  N  uptake 
between  the  two  cases  is  equivalent  to  the  total  amount 
ot  N  applied  in  one  wastewater  application.  These 
results,  therefore,  further  support  the  concept  that  inter¬ 
mittent  rainfall,  for  the  cases  studied,  directly  influences 
the  plant  uptake  rather  than  NO3-N  distribution  in  the 
soil  profile.  If,  however,  considerable  amounts  of 
NO3-N  are  present  in  the  upper  portion  of  the  soil  pro¬ 
file  prior  to  application  of  water  or  simulated  rainfall, 
profiles  of  NOj-N  vs  soil  depth  different  from  that 
shown  in  Figure  22  would  be  obtained. 
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DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

The  computer  program  that  we  developed  to  pre¬ 
dict  the  water  and  nitrogen  transport,  transformation, 
and  uptake  by  plants  in  die  soil  profile  in  land  treat¬ 
ment  systems  is  written  in  FORTRAN  language  and 
consists  of  a  SOURCE  (or  MAIN)  program,  10  sub¬ 
routine  programs,  three  FUNCTION  programs,  and  an 
input  data  section.  The  names  of  the  subroutines  are 
ID  1ST,  IDIST2,  ROOTS,  INDIDZ,  WATER,  WPROP, 
AMONIA,  NITRAT,  OUTPUT  and  TRIDM.  The  names 
of  the  function  programs  are  ZZ1 ,  ZZ2,  and  ZZ3.  In 
addition,  the  program  uses  subroutine  QSF  which  is  a 
standard  integration  subroutine  available  from  the  IBM 
System  360  Scientific  Subroutine  Package  (1970).  This 
subroutine  is  on  file  in  the  CRREL  computer  system. 

Input  data  section 

Unless  otherwise  stated,  the  format  for  the  input 
data  is  8F10.3  (see  FORMAT  100  in  the  SOURCE 
program).  The  input  data  which  must  be  provided  by 
the  user  are  as  follows: 

1 .  The  first  record  (card)  of  the  data  set  contains: 
DTT  =  initial  approximation  of  Ar 

DZZ  =  initial  approximation  of  (cm) 

2.  The  second  record  of  the  data  set  contains: 
SFLUX  =  flux  vof  wastewater  application 
(cm/h) 

ET  =  evapotranspiration  rate  T  (cin/h| 

QM  =  /max  of  the  Michaelis-Mentcn  equation 
for  nitrogen  uptake  (pg/h  cm  o(  root  length) 

QK  =  Michaclis  constant  Km  (pg/ml) 

CSNH4  =  concentration  of  NH4-N  in  the 
wastewater  Cs  (pg-N/ml) 

CSNOj  =  concentration  of  NO3-N  in  the  waste- 
water  Vs  (pg-N/ml) 

DISP  =  solute  dispersion  coefficient  D  (cm-'/hl 

3.  The  third  record  of  the  input  data  set  contains: 
CL  =  total  length  (cm)  of  the  soil  profile  L 
CL1  =  soil  depth  (cm)  to  the  first  soil  layer 

/  I  (see  Fig.  1) 

CL2  =  soil  depth  (cm)  to  the  second  soil  layer 
A ,  +/L2  (sec  Fig.  I). 

4.  The  fourth  record  contains  soil  water  parameters 
for  the  first  soil  layer  (see  Sensitivity  Analysis). 
The  format  used  here  is  E10.4,  6FI0.5  (sec 
FORMAT  500  in  SOURCE  program): 

AC1  =  tj  of  cq  37 
BCI  -  a  of  eq  37 
ATI  -  a  of  cq  36 
BT1  =  b  of  eq  36 

5.  This  record  is  similar  to  the  one  above  except 
for  the  second  soil  layer. 

6.  This  record  is  similar  to  the  one  above  except 
for  the  third  soil  layer. 


7.  This  record  contains  the  soil  bulk  density 
(g/cm3)  and  saturated  water  content  (cm3/cm3) 
for  each  soil  layer: 

ROU1  =  p|,  bulk  density  of  the  first  soil  layer 
THS1  =  (<>5)3 ,  saturated  water  content  of  the 
first  soil  layer 

ROU2  =  p 2,  bulk  density  of  the  second  soil 
layer 

THS2  =  (0S) 2,  saturated  water  content  of  the 
second  soil  layer 

ROU3  =  p 3,  bulk  density  of  the  third  soil  layer 
THS3  =  (0S) 3,  saturated  soil  water  content  of 
the  third  soil  layer. 

8.  This  record  contains  the  nitrogen  parameters 
for  the  first  soil  layers: 

RE XI  =  (X^)],  ammonium  distribution  co¬ 
efficient  (cm3/g)  for  the  first  soil  layer. 

RNIT1  =  (/e j  )| ,  nitrification  rate  coefficient 
(IT1 )  for  the  first  soil  layer. 

RDNIT1  =  (62)i,  denitrification  rate  coefficient 
(IT1 )  for  the  first  soil  layer. 

9.  This  record  is  similar  to  the  one  above  except 
for  the  second  soil  layer. 

10.  This  record  is  similar  to  the  one  above  except 
for  the  third  soil  layer. 

11.  This  record  contains: 

TINF  =  T,  duration  (h)  of  wastewater  (or  rain- 
f oil /  application 

TCYC  =  duration  of  wastewater  cycle  in  hours 
NCYC  =  number  of  cycles  desired. 

12.  This  record  contains: 

TW RITE  =  time  (111  at  whitb  output  data  are 
requested. 

Source  program 

The  main  functions  of  the  SOURCE  (or  main)  pro¬ 
gram  are  prescribing  the  DIMENSION  statements, 
reading  the  input  data,  and  setting  up  the  entire  sequence 
of  the  program.  A  detailed  flow  chart  of  the  SOURCE 
program  is  shown  in  Figure  24. 

The  DIMENSION  statements  arc  introduced  in  COM¬ 
MON  blocks  and  arc  labeled  LI  to  1.18.  The  array  si/c 
of  most  variables  is  set  to  be  3 1 0.  This  array  size  may 
be  changed  by  the  user  depending  on  the  depth  of  the 
soil  profile  ( L )  as  well  as  the  incremental  distances  A/. 

The  source  program  also  reads  the  input  parameters 
and  provides  the  output  for  these  parameters. 

The  following  variables  are  used  in  the  source  pro¬ 
gram: 

C  -  NH3-N  concentration  of  soil  solu¬ 
tion  C  (pg-N/ml),  dimension  =  NP1. 
Y  =  NO3-N  concentration  in  soil  solu¬ 
tion  Y  (pg-N/ml),  dimension  -  NP1. 
H  =  soil  water  pressure  head  h  (cm), 
dimension  =  NP1 , 
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Figure  24.  Detailed  flow  chart  of  simplified  nitrogen  model. 


TH  =  soil  water  content  0  (cm3/cm3), 
dimension  =  NP1. 

CON  =  soil  water  hydraulic  conductivity 
K(0)  (cm/h),  dimension  =  NP1 . 

CAP  =  soil  water  capacity  Cap  (h)  (cm'1 ), 
dimension  =  NP1. 

RDIST  =  root  density  distribution  R(z)  (cm), 
dimension  =  NX 

AA,  BB,  CC,  R  =  dummy  variables  which  are  used 
in  solving  the  matrix-vector  equa¬ 
tions  for  water,  NH4-N,  and  N03 
-N,  dimension  =  N 

XXX  =  soil  depths  (cm)  for  which  initial 
distributions  are  provided,  dimen¬ 
sion  =  NIN 

Cl  =  initial  distribution  of  pressure  head 
(cm)  dimension  =  NIN 
C2  =  initial  distribution  of  water  content 
(cm3/cm3)  dimension  =  NIN 
C3  =  initial  distribution  of  NH4-N, 
(/tg-N/ml)  dimension  =  NIN 
C4  =  initial  distribution  of  N03-N, 
(pg-N/ml),  dimension  =  NIN. 

Subroutine  IDIST 

This  subroutine  is  labeled  the  “Initial  Distribution 
Program”  in  the  program  listing  (Appendix  A).  This 
subroutine  uses  the  initial  input  distributions,  i.e.  initial 
conditions,  at  any  number  of  soil  depths.  In  order  to 
calculate  the  pressure  head,  water  content,  NH4-N  and 
NO3-N  concentration  distributions  at  incremental  dis¬ 
tances  DT  (CM)  for  the  entire  soil  profile  the  calcula¬ 
tions  are  carried  out  by  linear  interpolation  using  the 
input  data  points. 

In  this  subroutine  it  is  not  necessary  to  provide  initial 
water  content  distributions  corresponding  to  the  soil 
suction  for  the  various  soil  depths  (in  input  data).  Con¬ 
version  from  suction  or  pressure  head  to  water  content 
is  carried  out  in  subroutine  WPROP.  However,  the 
user  must  provide  some  value  (zero  is  recommended) 
for  all  input  data  points. 

It  should  be  emphasized  here  that,  if  input  data  for 
pressure  heads,  NH4-4,  etc.,  are  not  provided  for  the 
same  soil  locations  (depths),  the  user  must  use  subrou¬ 
tine  IDIST2. 

Subroutine  IDIST2 

This  subroutine  is  labeled  the  “Initial  Distribution 
Program  Number  2”  in  the  program  listing.  This  sub¬ 
routine  is  similar  to  subroutine  IDIST  except  that  it 
allows  the  calculation  of  initial  distribution  regardless 
of  locations  (soil  depths)  at  which  measurements  are 
provided.  Therefore,  this  subroutine  must  be  used  if 
initial  distributions  of  ail  variables  (pressure  head,  NH4 
-H,  etc.)  are  not  available  at  common  soil  depths. 


In  this  subroutine  each  variable  is  treated  separately. 
In  addition  the  number  of  data  points,  locations  and 
values  of  each  variable  must  be  supplied  in  the  main 
program. 

Subroutine  ROOTS 

This  subroutine  is  labeled  the  “  Root  Distribution 
Program"  in  the  program  listing.  The  subroutine  pro¬ 
vides  the  root  distribution  in  the  soil  profile.  This  root 
distribution  can  be  expressed  as  a  function  of  time  as 
well  as  soil  depth  in  the  profile.  This  is  provided  as  a 
mathematical  expression.  In  the  example  below,  an 
exponential  decay  function  of  root  length  with  soil 
depth  is  used  for  all  times.  If  desired,  the  user  can  ex¬ 
press  the  root  distribution  as  a  function  of  soil  depth  as 
well  as  time. 

Subroutine  INDTDZ 

This  subroutine  is  labeled  the  “Program  for  Adjusting 
Zone  of  Infiltration"  in  the  program  listing.  The  sub¬ 
routine  may  be  used  only  if  the  water  flux  in  the  soil 
profile  is  extremely  small.  In  such  a  case  it  is  reasonable 
to  be  mainly  concerned  with  the  top  portion  of  the 
soil  profile  during  the  initial  stages  of  wastewater  ap¬ 
plication  or  rainfall.  Such  an  assumption  is  applicable 
if  water  redistribution  continues  for  several  days  with 
no  new  wastewater  application  or  rainfall.  In  this  pro¬ 
gram  an  initial  length  of  30  cm  (root  zone)  is  assumed 
(N  =  30)  in  the  main  program.  This  length  is  automati¬ 
cally  increased  during  wastewater  application.  At  the 
termination  of  infiltration  the  total  length  of  the  pro¬ 
file  is  incorporated. 

It  is  important  to  emphasize  here  that  the  use  of 
such  a  subroutine  is  not  essential.  However,  its  use 
saves  j  considerable  amount  of  computer  time  especially 
during  simulation  of  infiltration  when  DT  is  smallest. 

If  this  feature  is  not  desirable,  the  user  may  ignore  it 
by  replacing  N  =  30  by  N  =  CL/DZ+0.01  in  the  main 
program  and  deleting  all  call  INDTDZ  statements. 

Subroutine  WATER 

This  subroutine  is  labeled  the  “Water  Flow  Pro¬ 
gram  ”  in  the  program  listing.  The  subroutine  provides 
the  solution  for  the  water  flow  equation  for  a  homo¬ 
geneous  or  a  layered  soil  profile.  The  method  of  solu¬ 
tion  is  based  on  the  finite  difference  approximations 
discussed  previously.  The  bottom  boundary  condi¬ 
tion  is  incorporated  in  the  solution  and  is  applicable 
for  a  water  table  boundary  or  a  semi-infinite  (i.e. 
deep)  soil  profile.  The  surface  boundary  condition 
(flux  type)  is  provided  from  main  program. 

Subroutine  WPROP 

This  subroutine  is  labeled  the  “Soil-Water  Properties 
Program”  in  the  program  listing.  The  subroutine 
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provides  the  soil-water  properties  for  each  soil  layer  in 
the  soil  profile,  namely  the  hydraulic  conductivity  as 
a  function  of  water  content  or  suction  and  the  water 
content-suction  relationship.  From  the  latter,  the 
water  capacity  term  is  calculated.  In  this  example  (see 
Sensitivity  Analysis)  mathematical  expressions  are  used 
to  describe  these  relationships.  This  subroutine  is 
called  by  subroutine  WATER  for  every  time  step. 

Subroutine  AMONIA 

This  subroutine  is  labeled  the  “Ammonium  Trans¬ 
port  and  Transformation  Program”  in  the  program  listing. 
The  subroutine  provides  the  solution  to  the  ammonium 
transport  and  transformation  equation  under  transient 
flow  conditions.  It  also  calculates  the  ammonium 
uptake  by  plant  roots.  The  method  of  solution  is  the 
finite  difference  approximation  method.  The  rate  of 
nitrification  and  denitrification  and  the  distribution 
coefficient  for  NH4-N  release  are  obtained  from  func¬ 
tions  ZZ1,  ZZ2,  and  ZZ3,  respectively  (see  below). 

Subroutine  NITRAT 

This  subroutine  is  labeled  the  “Nitrate  Transport  of 
Transformation  Program”  in  the  program  listing.  The 
subroutine  provides  the  solution  to  the  nitrate  trans¬ 
port  and  transformation  equation  under  transient  flow 
conditions,  (t  also  calculates  the  nitrate  uptake  by 
plant  roots.  The  rate  of  nitrification  and  denitrifica¬ 
tion  and  the  distribution  coefficient  for  NH4-N  ex¬ 
change  are  obtained  from  functions  ZZ1 ,  ZZ2,  and 
ZZ3,  respectively. 

Functions  ZZ1,  ZZ2,  ZZ3 

Function  ZZ1  provides  the  rate  coefficient  for  nitri¬ 
fication  (£j )  for  each  soil  layer.  In  the  example  de¬ 
scribed  (see  Sensitivity  Analysis)  the  rate  coefficient 
was  considered  dependent  on  the  soil  water  suction. 

It  was  assumed  that  nitrifying  bacteria  are  present  in 
optimum  number  and  the  change  in  their  population 
during  a  cycle  of  wastewater  application  (most  often 
weekly)  is  negligible. 

Function  ZZ2  provides  the  retardation  factor  R  for 
ammonium  exchange.  This  is  achieved,  for  each  soil 
layer,  from  KD,  0  and  p  at  incremental  soil  depths. 

All  the  above  functions  (ZZ1,  ZZ2,  and  ZZ3)  are 
called  by  subroutine  AMONIA  and  NITRAT  at  every 
time  step  DT. 

Subroutine  OUTPUT 

The  main  function  of  ibis  subroutine  is  to  print  the 
results  (output  data)  at  specified  times.  A  second  func¬ 
tion  of  this  subroutine  is  to  carry  out  several  integration 
in  order  to  calculate  the  total  amount  of  N03-N  and 
NH4-N  in  the  soil  solution  and  total  NH4-N  in  the  ex 
changeable  phase. 


Subroutine  TRIDM 

This  subroutine  is  labeled  the  “Tridiagonal  Matrix 
Program”  in  the  program  listing.  The  subroutine  pro¬ 
vides  the  solution  of  a  tridiagonal  matrix-vector  equa¬ 
tion  (see  text).  This  subroutine  is  utilized  by  subroutine 
WATER,  AMONIA,  and  NITRAT  at  every  time  step 
DT. 
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APPENDIX  A:  PROGRAM  LISTING 


C*** **************************************************** »*»**••*****••*****»*•*£ 

c  c 

C  A  SIMPLIFIED  MODEL  FOR  PREDICTION  C 

C  C 

C  OF  NITROGEN  BEHAVIOR  IN  LAND  C 

C  C 

C  TREATMENT  OF  BASIE  NATER  C 

C  C 

c** *****«***«*•*************»***•»****•»**•*»•«»••*•*•*#»*»*•••» **********»****c 

c  c 

C  TdE  PURPOSE  OF  THIS  MODEL  IS  TO  PREDICT  TUB  BEHAVIOR  OF  THE 

C  AMMONIUM  AND  NITRATE  NITROGEN  SPECIES  IN  THE  SOIL  IN  LAND  TREATMENT  C 

C  SYSTEMS..  THE  PROGRAM  IS  UASED  ON  THE  SOLUTION  OF  THE  TRANSIENT  C 

C  WATER  FLOW  EQUATION  SIMULTANEOUSLY  WITH  THE  EQUATIONS  DESCRIBING  C 

C  THE  TRANSPORT  ,  TRANSFORMATIONS,  AND  UPTAKE  OF  NITROGEN  BY  PLANTS.  C 

C  C 

C  MAIN  FEATURES:  C 

C  THE  MODEL  IS  VALID  FOR  UNIFORM  AS  WELL  AS  MULTILAYERED  SOIL  C 

C  PROFILES.  IN  ADDITION,  ThE  PROGRAM  IS  FLEXIBLE  AND  IS  DESIGNED  C 

C  TO  INCORPORATE  THE  FOLLOWING  (INPUT)  CONDITIONS  AS  DESIRED  ;  C 

C  1.  RATE  OF  WASTE  WATt.fi  APPLICATION  C 

C  2.  DURATION  CF  WASTE  WATER  APPLICATION  C 

C  3.  DtPTu  OF  INDIVIDUAL  SOIL  LATtfiS  C 

C  4.  CONCENTRATION  OF  AMMONIUM  AND  NITRATE  iN  THE  C 

C  WASTE  WATER  C 

C  5.  WASTE  WATER  APPLICATION  CYCLE,  l.E.  SCHEDULING  C 

C  6.  SOIL  WATER  PfiOPEFTIES  AND  NITROGEN  TRANSFORMATION  C 

C  MECHANISMS  FOR  INDIVIDUAL  SOIL  LAYERS  C 

L  7.  PLANT  ROOT  DlSIfi  I3UTI0N  AND  JKO'WTu  IN  IHE  SOIL  C 

C  H.  HATE  CF  NITROGEN  UPTAKE  BY  PLANTS  C 

C  9.  FVAPOTSANSPIR ALION  RATE 

C  10.  INITIAL  DISTRIBUTION  OF  WATER  ANC  NITROGEN 


C  SPECIES  IN  ThE  SOIL  PROFILE 

C WWW********************************************************************** ***w*c 

c 

COMMON/ LI/  C  (310) , Y  (31C) 

COMMON/L2/  A A (310)  , 30(310) ,CC (310) ,R (310) ,KDIST (31 C) 

COMMQN/L3/  N,NM1,NK2,NF1,NP2 
COMMGN/L4/  ALPHA, BETA, DT,DZ 
COKMCN/L5/  NX,NX1, N?MAX,CCN1 
C0.1MON/L6/  FLNH4,  PLNO.I,  DNIIRF 

COMMON/L7/  SFLU  X,  ET,QM,  QK,CSNIi«,  C5N03,  DISP,XL 

COMMCN/Lh/  XXX (30) ,C1 (30) ,C 2 (30)  ,C3  (30)  ,C4  (30)  ,NIN 
COMMCN/L9/  TIME, TINE, 7CYC 

COMM JN/L  1 0/  H  (  3  1 0)  ,  CON  (310)  ,CAP(310)  ,TM(310) 

COMMON/L 1 1/  CL,CL1,CL2,L1,L2 
COMMON/L  1  2/  AC  1 , DC  1  ,  AT  1 , B1 1 
CO MMC  N/l.  1 3/  AC2,BC2,  AT2,bT2 
COMMO N/L  1 4/  AC3,EC3,AI3,2T3 

COMMON/L 1 5/  ROU 1 ,THS 1 , ROU2, THS2 ,ROU? ,THS3 
C0MM0N/L10/  REXl,RNIT  1,RCNIT1 
COMMCN/L17/  ELX2,RN1T2, RLNIT2 
CO MMON/L 10/  R£X  3, P.NIT3,  RDNIT 3 
ICO  FOPMAT  (HF1J. 3) 

101  FORMAT (5X, 'INITIAL  DT,  Hu  = • , FI 0 . S/, 5X , • INI II AL  DZ  , 

200  FORMAT  (2F1C.5,13) 

300  FORMAT  (214) 


cM=',no.  V/) 


n  n 


600  FORMAT  (10E12.4) 

500  FORMAT  (E10.4,6P10.5) 

105  FOBBAT  (SX.'FLUX  OF  HASTE  HATEB  APPLICATION,  CB/  HB  =',P10.5/, 

*5X, ' E VAPOTBANSPIBATION  FATA,  CM/HR  =',F10.5/, 

I5X, 'NITROGEN  UPTAKE  BATE,  BICROGBAB-H/CB  OF  BOOT  LENGTH  PEE  HOUR', 
$•  =  •  , F 10. 5/, 

*5X, 'MICHAELIS  CONSTANT  =',F10.5/, 

$5X,  •CONCENTRATION  OP  APPLIED  NH4-N  ,  HG/LITHE  =*,F10.5/, 

*5X, 'CONCENTRATION  OF  APPLIED  N03-N  ,  BG/LITBE  =',P10.5/, 

I5X, 'SOLUTE  DISPERSION  COEFFICIENT,  CN**2/HH  =',F10.5/) 

110  FOBBAT  (5X, 'TOTAL  LENGTH  Ot  SOIL  PBOFILE,  CB  =',P10.5/, 

♦5X , • SOIL  DEPTH  TC  THE  FIRST  SOIL  LAYER,  CB  =',F10.5/, 

J5X , ' SOIL  DEPTH  TO  THE  SECOND  SOIL  LAYER,  CB  =',F10.5/) 

115  FOBBAT  (5X, 'SOIL  NATES  PAHAHETERS  FOB  THE  FIBST  LAYER  :',4E15.4) 

120  FOBBAT  (5X, 'SOIL  HATES  PARABETERS  FOR  THE  SECOND  LAYER  :',4E15.4) 

125  FORMAT  (5X,'SCIL  HATER  PARABETERS  FOR  THE  ThIRD  LAYER  :',4E15.4) 

130  FORMAT  (5X, 

♦'FIRST  LAYER  ;  BULK  DENSITY  =• ,P 1 0. 5, 10X , 'SATURATION  =',F10.5/, 
♦5X, 

♦•SECOND  LAYER  ;  BULK  DENSITY  =• ,F 1 0. 5 , 1  OX ,' SATURATION  =',F10.5/, 
*5X, 

♦•THIRD  LAYER  ;  BULK  DENSITY  =', F 10. 5, 10X, • SATURATION  =',F10.5/, 

♦  5  X ) 

135  FORMAT  (5X, 'FIRST  LA Y ER : • , 1 0X, • Hh4- N  EXCHANGEABLE  COEFFICIENT,  CB3/ 
♦GB  =',F10. 5/, 26X, 'NITRIFICATION  RAT F  COEF., HR-1  = •  , FI 0. 5/.26X, 
♦•DENITRIFICATION  RATE  COEF.,  HR- 1  =',E10.5) 

140  FORBAT  (5X, 'SECOND  L A YER • , 1 0X , » NH4- U  EXCHANGEABLE  COEFFICIENT,  CB3/ 

♦  GB  =',F10. 5/, 26X, 'NITRIFICATION  RATE  COEF. , HR- 1  =• , F  1 0. 5/, 26X, 
♦•DENITRIFICATION  RAT F  COEF. ,  HR- 1  =',F10.5) 

145  FORMAT  (5X, 'TRIED  LA Y ER : • , 1 0X , • NH4- N  EXCHANGEABLE  COEFFICIENT,  CB3/ 

♦  GB  =',F1 J.S/,26X, 'NITRIF ICATION  RATE  COEF., HR-1  *•  ,  f  10. 5/, 26X, 
♦'DENITRIFICATION  RATE  COEF.,  HR-1  =',F10.5) 

150  FORBAT  (5X, 'DURATION  OF  HASTE  HATER  APPLICATION  ,  HRS  = • , F 10. 5/, 5X  , 
♦'SCHEDULE  OF  HASTE  HATER  APPLICATION,  I.E.  CYCLE  DURATION  =•  , 
tFIO. 5/, 5X, 'NUMBER  OF  CYCLES  =',I3/) 

155  FORMAT  (5X, 'TIBE  AT  WHICH  OUTPUT  DATA  IS  REQUESTED, HE  =',F10.5/) 

160  FORBAT  (125, ' I  N  P  U  T  DAT  A*,///) 

165  FORBAT  ( '  1 ' ) 

WRITE  (6,160) 

READ  (5,100) 

WRITE (6,101) 

READ  (5, 100) 

WRITE  (6,  105) 

READ  (5,100) 

WPITE  (6, 110) 

READ (5,500) 

WPITE  (n, 115) 

READ  (5,500) 

WRITE  (o,  120) 

READ  (5,500) 

WRITS  (6,  125) 

READ  (5, 100) 

•  RITE  (b, 130) 

RF AD  (5,  100) 

WRITS  <6,  1  J'j) 

READ  (5,  100) 

•  RITE  (fc, 140) 

R  EAD  (5,  100) 


DTT,DZZ 
DTT, DZ? 

SFLUX, FT,vB,CK»CSNH«,CSNOJ,DISP 

SFLUX, ST.^B, 0K,CSNH4,CSSO3, DISP 

CL,CL 1 ,CL2 

CL, CL  1 ,CL2 

AC  1 ,  DC  1 ,  AT  1 ,  BX  1 

AC  1 ,  BC 1  ,  AT  1 ,  13T1 

AC2,PC2,AI'2,BT2 

AC2,oC2, AT2,bT2 

ACJ,BC3,AT3,ET3 

AC3.BC3, AT3,BT3 

FOU1,  TtiS  1,bO'J2,THS2,ROU3,iHS3 

ROU1,THS1,POU2,THS2,ROU3,TH33 

FEX1,RNIT1,RCNIX1 

REXl,RNIT1,IiDNIT1 

REX2,SNIT2,RDNIT2 

REX2,r<NII2,RCNIT2 

REX3,RNI13,RD.UT3 
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c 


c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


l. 


VBITE  (6,  145) 
HEAD  (5,200) 
WRITE  (6,  150) 
BEAD  (5,100) 
SHITE  (6, 155) 


BEX3,HNIT3,RDNIT3 

1INF,TCYC,NCYC 

TINF,TCTC, NCYC 

1BBITE 

TSBITE 


FLN03*0. 0 
PLNH4=0.0 
DNITRF-0.0 


DT=DTT 

DZ=DZZ 

fcFLX=SFLUX 


NKK=TWBITE*6,  10 
N=CI/DZ*0. 10 
L1=CL1/DZ*0.  10 
L2=CL2/DZ+0.  10 
Ntt1  =  N-  1 
NH2=N-2 
NP1=Nt  1 
BP2=N*2 


ICL=CI/DZ+0.  10 
1CL-ICL*  1 


BEAD  NUHBEB  (INTEGER)  OF  DATA  POINTS  FOB  BtilCH  INITIAL  DISTRIBUTIONS 
APE  PROVIDED 
BEAD  (5, 330)  N1N 

BEAD  SOIL  DEPTHS  FOR  INITIAL  DISTRIBUTION 
BEAD  (5,  1  00)  (XXX(I) ,I=1,NIN) 

BEAD  PiESSIBE  HEAD  FOB  INITIAL  DISTRIBUTION 
BEAD  (5,  100)  (Cl  (I)  ,1=1, NIN) 

PEAO  JATEH  CONTENT  FOB  INITIAL  DISTBI3UT ION  (IF  AVAILABLE) 

READ(5,100)  (C2  (I) ,1=1, NIN) 

BEAD  NB4-N  COMCEN,  FOB  INITIAL  DISTRIBUTION 
BFAD  (5,  1  OC)  (C3  (I)  ,  1  =  1  ,  N  I N) 

READ  N03-N  CONCEN,  FOB  INITIAL  DISIR1CUTION 
H  E  A  D  (  5 ,  100)  (C4  (I)  , 1=1 , N IS) 

CALL  IDIST 


IF  A-L  VARIABLES  (  l.E.  PRESSURE  HEAD,  NK4-N  CONC.  ETC.)  ABE  NOT 
AVILAELE  AT  CCflfiCN  SOIL  DEPTHS,  SUBROUTINE  IDIST2  1UST  BE  USED. 
HEAD  NUMBER  OF  DATA  POINTS  FOR  PRESSURE  HEADS 
BEAL  (5, JOC)  SIN 

READ  SOIL  DFt-TH  (LOCATIONS)  FOB  PRESSUBE  HEADS 
BEAD (5, 1 OC)  (XXX (I) ,I=1,NIN) 

BEAD  PRESSURE  HEAD  FOB  COER ESPON DIN 0  DEPTHS 
R I  A  D  (  5 , 100)  (Cl  (I)  ,1*1, NIK) 

CALL  I D 1ST  2 
DO  11  1=1, NT1 

11  H  (I)  =R  (I) 

LEAD  NUHBEB  OF  DATA  POINTS  FOB  WATER  CONTENT  (IF  AVAILABLE) 

READ  (5, JOQ)  NIN 
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n 


C  BEAD  SOIL  DEPTH  (LOCATIONS)  FOB  WAT  Eli  CONTENT 

BEAD  (5, 100)  (XXX  (I) ,1=1, NIN) 

C  BEAD  NATEB  CONTENT  FOB  CORRESPONDING  DEPTHS 

BEAD  (5, 100)  (C 1  (I) ,  1  =  1 ,  NIN) 

CALL  ID1ST2 
DO  12  1*1, NP1 

12  TH  (I)  =R  (I) 

C  BEAD  N'JSBER  OF  DATA  POINTS  FOB  NH4-N  CONCENTRATION 

BEAD  (5, 300)  NIN 

C  BEAD  SOIL  DEPTH  (LOCATIONS)  FOB  NH4-N  CONCENTRATION 

READ  (5,  100)  (XXX  (I)  ,1  =  1, NIN) 

C  BEAD  NhU-N  CONCiN.  FOB  CORRESPONDING  DEPTHS 

BE AD  (5 , 1 30)  (Cl  (I)  ,1=1, NIN) 

CALL  IDIST2 
DO  13  1=1, NF1 

13  C  (I)  =  E  (I) 

C  BEAD  NUM3ER  GF  DATA  POINTS  FOR  N03-N  CONCENTRATION 

BEAD  (5, 30C)  NIN 

C  BEAD  SOIL  DEPTH  (LOCAXIONS)  FOB  N03-N  CONCfcNTBATIO N 

READ  (5,  100)  (XXX  (I) ,1=1, NIN) 

C  READ  N03-N  CCNCSN.  FOB  CORRESPONDING  DEPTHS 

READ  (5, IOC)  (Cl  (1) ,1=1, NIN) 

CALL  IDIS12 
DO  14  I =  1 , N  F 1 
1  4  Y  (I)  =8  (I) 


C 


C 


c 


5 


C 


6 


CALL  BOOTS 
CALL  WPfiOP 
CALL  CUTFUT 
WRIT  E  (€,165) 

T1SE=0. 0 
XTINF=TINF 

DO  39' j  IJKL=1,NCYC 

TINF=TJME*X7INF 

DT=DTT 

DZ=DZZ 

SFLUX=WFIX 


N  =  30 
NX  =  N 

CALL  INDTD2 
NX1=NX-1 


IF  (SFLUX.L2. 0.0)  GO  TO  26 
TTT  =  4. 0/  (SFLUX* 1000.  0) 

IF  (DT.LK.T7T)  GO  TO  5 
NT=DT/TTI»0,  10 
DT  =D i/NT 
CCJTINUK 

AIP.IA  =  DT/  (2. 0*D2*DZ) 

LE?A=D1/DZ 

CALL  WIFfJP 

IF  (CUNl.Ge.  (SFL’JX/2))  GO  TO  15 
6  (1)  =:i  (1)  <-20 
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CALL  WPKCP 

IF  (C0M1.lt.  (SFL0X/2) )  GO  TO  6 
1F(H (1) .GT.O.O)  H  (1) -0. 0 
H (1) =H  (1 > -60. 0 
C 

XL=  1000.0 
DELli  =  2 . 0 
7  CONTINUE 
IL2  =  30 

DO  25  IK=  1 ,  ILZ 
CALL  MATES 
CALL  AMONIA 
CALL  NITDAT 
H (1) =a  (1) +  DELH 
25  CONT I N  US 

TI'1E  =  Ti;i£*DT*ILL 
CALL  INDTD2 
C 

15  CONTINUE 
C 

DT  =  QT  *2 . 0 

ALPHA=CT/  {2. 0*DZ*DZ) 
aFTA=ET/D7 
XI=TINF-riHE 
XT  =  Xl/5 

IL=XT/DT»0.010 
JO  30  IZFT=1,5 
DO  1 7  1=1, NX 

17  AA  (I) =1DIST  (I) *CCN  (I) 

CALL  gSt’(DZ,AA,R,NX) 

XL=  3  (NX) 

DO  1C  11=  1, IL 

ADJ= JZ*  ( 1 . 0— SFLU X/CON 1) 

H (1) =H (2) -ADJ 
CALL  NATES 
CALL  A3CNJ.A 
CALL  NIT 3 AT 
10  CONTINUE 

TI.T2  =  Tl.-iE*IL*DT 
CALL  INDTD2 
30  CONTINUE 

K  K=T I M  E ♦ 0 • 0  5  C 
J  =  M\-NKE*  (KK/NKK) 

IF(J.tw.O)  CALL  OUTPUT 
SFLU  X  =  0. 0 

TN=6. 0+TCYC*  (1JKL-1) 

53  IF  (  (TN-TINF) -GT. 1.0)  GO  TO  54 

TS=T  N  +  fi.  0 
GO  TC  53 

54  iL=  (T N-T 1 1 E) /DT ♦ 0. 0  1  0 
JT=  (TN-T I TE) /I l 


DO  31  1=1, NX 

31  AA (I) =LD1ST  (I) *CON  (1) 
CALL  CSF  (D7. ,  A  A  ,  ft ,  NX) 

X  L  =  n  (NX) 

DO  5  2  il=  1 ,  IL 
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ADJ  =  DZ*  ( 1 . Q-SPLUX/CON 1 ) 

H  ( 1) ®H (2) -ADJ 
CALL  WATER 
CALL  AHCNIA 
CALL  N1TNAT 
52  CONTINUE 

TIHE=TINE*tL*DT 
CALL  INCTDZ 
KK=TIHE*0.050 
J=KK-NKK*  (KK/NKK) 

If  (J. EQ. 0)  CALL  OUTPUT 
C 
C 

26  CONTINUE 
C 

DT=DT*2 

IF (DT.JT.0.05G.AND.DT.LT.0.  10)  DT  =  0. 10 
ALPHA=UT/(2.0*DZ*DZ) 

Df TA=CT/D2 
IL=6.0C/DT*0.010 
DT  =  6. C/I L 
C 

DC  34  1 = 1 , NX 

34  AA (I) =RD1ST (I)*CON  (I) 

CALL  CSF (DZ , A  A,  R , NX) 

XL  =  I?  (NX) 

DC  36  11=1,11 
AEJ  =  CZ*  ( 1. 0-SFLUX/CON 1) 
a  11)  *0  (2)  -  ADJ 
CALL  KATtfl 
CALL  AHCNIA 
CALL  NITNAI 
3b  CONTINUE 

TIHE=TIME»IL*DT 

KK=TIHE*0.050 

J  =  KK- NKK*  (KK/NKK) 

IF(J.EQ.C)  CALL  CUTPUT 

DT  =  2*DT 

IF  (DT.UT. 0.050. AND. D'i.LT.O.  10)  DT=0.10 

IF  (DT.OT.C. 100. AND.DT.LT. 0.40)  DT  =  0.250 

AL?HA  =  0T/  (2.0*DZ*CZ) 

L'EiA  =  Cl/DZ 

TIHH=T1ME-TCYC*  (1JKL-1) 

17.  EFT  =  (TCYC-TIHH)/6.C 
DO  40  1ZZ=1,IZEFT 
IL=6. 00/DT+0.010 
JT  =b . C/IL 
DO  18  1=1, NX 

10  A  A  (I)  =!.D  1ST  (I )  *CON  (I) 

CALL  CSF  (DZ, AA,b,NX) 

XL  =  j  (NX) 

DO  35  11=1, XL 

AD. 1  =  02*  (1.0-KFLU X/CON1) 

li  (1)  =.l  (2)  -ADJ 

CALL  WATER 
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4 


f 


1 


CALL  AMONIA 
CALL  N1TEAT 
35  CONTINUE 

T1ME=T1HE*IL*DT 
KK=TIHE*0. 050 
J  =  KK-NKK*  (KK/NKK) 

If  (J.  EQ.  0)  CALL  OUTPUT 
DT=2*DT 

IF  (DT •  GT. 2. 0)  DI  =  2.  0 
ALPHABET/  (2*EZ*DZ) 

BET  A=  CT/DZ 
40  CONTINUE 

C  IF  ANY  CHANGES  IN  INPUT  DATA  (  ESPECIALLY  THE  BOUNDARY  CONDITIONS) 

C  ARE  REQUIRED  IN  THE  CYCLE,  THE  INPUTS  SHOULD  BE  ENTERED  HERR 

C  FOB  EXAMPLE  NEK  FLUX,  EVAPOT RANSPIkATION  RATE,  DURATION  OF  HASTE 

C  HATER  APPLICATION,  CYCLE  DURATION,  CONCENTRATION  OF  NH4-N ,  AND 

C  NO  3-N  MAY  3E  ENTERED  AS  FOLLOWS. 

C  READ  (5, 100)  SFLUX,  ET , T I N F , TC YC , CS N H4 , CS NO i 

C  t.S.  RAINFALL  EVENTS  CAN  BE  TREATED  AS  A  CYCLE  BY  ENTERING  ZERO 

C  VALUES  FOE  APPLIED  NH4-N  AND  N03-N  AS  HELL  AS  PROVIDING  THE  PROPER 

C  INPUTS  FOR  FLUX  ( I NT ENSIT Y) , DUR AT10N  OF  RAINFALL  AS  HELL  AS  THE 

C  TOTAL  TIME  BEFORE  THE  NEW  WASTE  WATER  APPLICATION  BEGINS  (TCIC) 

C  WARNING: 

C  IF  THE  ABOVE  IS  DESIRABLE  THEN  NEW  DATA  FOR  EACH 

C  SUBSEQUENT  CYCLE  MUST  BE  PROVIDED 

J  9  9  CONTINUE 
C 

STOP 

END 
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£**•*•*****»•**••••*•*»**•«********••***•»•**•*•»••••*•*•*••*•***•*•***••< 

c 

C  I  N  1  I  A  L  DISTRIBUTION  PROGRAM 

t 

C  THIS  SUDRCj  JI INE  UTILIZES  THE  INITIAL  INPUT  DISTRIBUTIONS  ,1.E. 

C  INITIAL  CONDITIONS,  AT  AN  *  NUMBER  OF  SOIL  DEPTHS  IN  THE  SOIL, 

C  IN  ORDER  TO  CALCULATE  THE  PRESSURE  HEAD,  HATER  CONTENT,  NH4-N 

C  AND  N03-N  CONCENTRATION  CISTR ICUBIONS  AT  INCREMENTAL  DISTANCES  DZ 

C  (Cl)  FOR  THE  ENTIRE  SOIL  PROFILE.  IN  THIS  PROGRAM 

C  THE  CALCULATIONS  ARE  CARRIED  OUT  USING  LINEAR  INTEREOLATION 

C  USING  ThE  INFUT  DATA  FOINTS. 

C  IMPORTANT: 

C  IN  THIS  PROGRAM  IT  IS  NOT  NECESSARY  TO  PROVIDE  INITIAL 

C  RATER  CONTENT  DISTRIBUTIONS  CORRESPONDING  TO  TdE  SOIL 

C  SUCTION  FOR  THE  VARIOUS  SC1L  DEPTH  (IN  INPUT  DATA  ). 

C  CONVERSION  FROM  SUCTION  OR  PRESSURE  HEAD  TO  WATER  CONTENT 

C  IS  CARRIED  CUT  IN  SUBROUTINE  WPROP  HOWEVER,  THE  USER 

C  MUST  PHIVIDE  SOME  VALUE(ZERO  IS  RECOMMENDED)  FOB  ALL 

C  INFUT  DATA  POINTS 

C 

C  WARNING 

C  1?  INPUT  DATA  FOR  PRESSURE  HEADS,  NH4-4,  ETC.  ARE  NOT  PROVIDED 

C  FOR  TllE  SAME  SOIL  LOCATIONS  (DFPTIiS)  THE  USER  MUST  USE 

C  SUBROUTINE  1DI3T2 

C 

C***« ******************************** ***«**•******•***•****«*•*•**•**••*•*•*•**(; 

C  C 

SU3KGUT1N2  ICiST 
COMMON/ LI/  C  (310) ,Y  (31 C) 

CCMMON/L3/  K,NM1,NH2,NP1,NP2 
C011CS /LA/  ALPHA, fiETA,CT,D2 

COMMCN/LH/  XXX  (30)  ,C1  (30) ,C2(30)  ,C3  (30) ,C4(30) ,NIN 
COM  MON/E  10/  ti  (310)  ,CCN  (310)  ,CAP  (310)  ,TH  (310) 

COMHON/lll/  CL,CL1  ,CI.2,L1,L2 
1)3  FORMAT  (HE  1  3. 4) 

1=  1 

CO  2  C  K  = 1 , N  P  1 
A-JZ*  |K- V) 

5  IF  (A.IF.XXX  ( X ♦ 1 > )  GO  TU  10 
1  =  1*1 
ij C  TO  » 

10  ti(X)=Cl(l)*(A-XXX(I))*((C  1(1*1)  -Cl  (I)  )  /  (XXX  (1*1)  -XXX  (I)  )  ) 

1H (K) =C2  (I) ♦ (A-XXX (1) ) * ( (C2 (1*1) -C2 (I) )/  (XXX  (1* 1)-XXX  (I) ) ) 

C  (X)  =C3  (I)  ♦  (A-XXX  (I)  )  *  ( (CM  (1*1)  -C3  (I) )/  (XXX  (1*1)  -XXX  (I)  )  ) 

Y  (K)  =C  »  (I)  ♦  (A-XXX  (1)  )*  (  (C4  (1*1) -C4  (I)  )/  (XXX  (1*1)  -XXX  (I)  )  ) 

20  CONTINUE 
RETURN 
EN  J 
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nno  n  n  n  n  n  n  n  n  n  n  n  n  n  n  n  n  n  nn 


or  no  o 


r ***************************************************** **♦♦♦**♦♦**• *************0 

C  c 

C  INTIAL  DISTRIBUTION  PROGRAM  NO.  2  C 

C********** ************  **************  *****************  ******** *****************c 

C  14IS  SUBROUTINE  IS  SIMILAR  TO  SUBROUTINE  IOIST  FXCE’PT  THAT  c 

C  IT  ALLOWS  THE  CALCULATION  OF  INITIAL  DISTRIBUTION  REGARDLESS  OF  C 

C  LOCATIONS  (SOIL  DEPTHS)  AT  WHICH  MEASUREMENTS  ARE  FBOVIDED.  C 

C  THEREFORE,  THIS  SUBROUTINE  MUST  BE  USED  IF  INITIAL  DISTRIBUTIONS  C 

C  OF  ALL  VARIABLES  (  PRESSURE  HEAD,  NH4-N,  ETC.)  ARE  NOT  AVAILABLE  AT  C 

C  COMMON  SOIL  DEPTHS.  ~ 

C  IPIFOWTANTj  ^ 

c  EACU  VARIABLE  IS  TREATED  SEPENATLY.  IN  ADDITION  THE  C 

C  NUMBER  OF  DATA  POINTS,  LOCATIONS  AND  VALUES  OF  EACH  C 

C  VARIABLE  MUST  BE  SUPPLIED  IN  THE  MAIN  PROGBAM.  C 

^** ******** ******************* ***************** ******* *************************C 


SUBROUTINE  IDIST2 

COMNON/L  2/  A  A  (3  10)  , B3  (31 0) ,CC  (3 10) ,  R  (3 10) ,RDIST (31 0) 
COMHON/Li/  N, NH1, NM2, NE 1, NF2 
COMMON/LV  ALPHA, BETA, Cl, DZ 

COMMCN/Ld/  XXX  (30) ,C1  (30) ,C2  (30)  ,C3  (30) ,C4  (30)  ,NIN 
1=1 

DO  20  K  =  1 , N  F 1 
A  =  DZ*  (K-1) 

5  IF  (A. LE. XXX  (1*1) )  GO  TO  10 
1  =  1*1 

GO  TO  5 

10  I<  (K)  =C1  (1)  *  (A-XXX  (I)  )  *  (  (Cl  (1*1)  -Cl  (I)  )  /  (XXX  (1*1)  -XXX  (I)  )  ) 
20  CONTINUE 
RETURN 
END 


C  C 

C  ROOT  DISTRIBUTION  PROGRAM  C 

c  c 

C  THIS  SUBROUTINE  PROVIDES  THE  HOOT  DISTRIBUTION  IN  THE  SOIL  PROFILE.  C 

-  THIS  EOOI  DISTRIBUTION  CAN  BE  EXPRESSED  AS  A  FUNCICN  OF  TIM*1  AS  C 

C  WELL  AS  SOIL  DEPTH  IN  THE  PROFILE.  C 

0  THIS  MAY  MAY  PROVIDED  AS  A  MATHEMATICAL  EXPRESSION.  C 

C  IN  THE  EXAMPLE  BELOW  AN  EXPONENTIAL  DECAY  FUNCTION  CF  ROOT  LENGTH  C 

C  WITH  SOIL  DEPTH  IS  USED  EOF  ALL  TIMES.  C 

C  IMPORTANT :  C 

FOCI  DISTRIBUTION  IS  EXPRESSED  IN  TERMS  CF  TOTAL  ROOT  C 

LENGTH  PER  UNIT  BULK  VOLUME  CF  SOIL.  C 


SUBROUTINE  ROOTS 

COM  .MG  N/L2/  A  A  (310)  ,tE  (310)  ,CC  (310)  ,3  (310)  ,RDIST  (310) 

COMMON/ L  i/  X,  NM.Nf"2,NP  1,NF2 
COKMCN/L4/  ALPHA, PET  A, CT,CZ 
COMMON/LV  NX  ,  K  X 1 ,  NRMAX,CON1 

cony.cn/vt/  time.tinf  ,tcyc 

MATHEMATICAL  EXPRESSION  CF  ROOT  LENGTH  AS  A  FUNCTION 
Of  SOIL  DEPTH 
DO  5  I  =  1 ,  N 
Z~  (1-1)  *Ds 

b  «DIS T (I) =226. C* EXP  (-0. 1 0 0 * 2 ) 

C 

MX-  B  i) 

10  CT  =  PCIST  («X)/nP.IST  (1) 

IF (CT.LT. 0.010)  GO  1C  1b 

n  x = :j  x  ♦  1  o 

GO  TO  10 
15  NX  1  =  NX-  1 
VV.  '  A  X  =  N  X 
RETURN 
END 
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Q* ************************************************************************** ***Q 

c 

C  P  R  C  G  R  A  M  FOR  ADJUSTING  ZONE 

C 

C  OF  INFILTRATION 

C 

Q*********************  **************************************************** *****Q 


C 

C  THIS  SUBROUTINE  MAY  BE  USED  ONLY  IF  THE  HATER  FLU*  IN  THE  SOIL  C 
C  PROFILE  IS  EXTREMELY  SMALL.  IN  SUCH  A  CASE  IT  IS  REASONABLE  TO  C 
C  oE  MAINLY  CONCERNED  WITH  THE  TOP  PCRTCON  OF  THE  SOIL  PROFILE  DURING  C 
C  THE  INITIAL  STATES  OF  HASTE  HATER  APPLICATION  OH  RAINFALL.  C 
C  SUCH  AN  ASSUMPTION  IS  APPLICABLE  IF  HATER  REDISTRIBUTION  CONTINUES  C 
C  FOR  SEVERAL  DAYS  HIT .1  NO  NEW  HASTE  HATER  APPLICATION  OR  RAINFALL.  C 
C  IN  THIS  PROG  RAM  AN  INITIAL  LENGTH  OF  30  CM  IS  ASSUMED  (N  =  30).  C 
C  THIS  LENGTH  IS  AUTOMATICALLY  INCREASED  DURING  HASTE  HATER  APPLICATION.  C 
C  AT  THE  TERMINATION  OF  1NFIT  RAT  ION  THE  TOTAL  LENGTH  OF  THE  PROFILE  IS  C 
C  INCORPORATED.  C 
C  IMPORTANT:  C 
C  THE  USE  CF  SUCH  A  SUBROUTINE  IS  NOT  NESESSARY.  C 
C  HOWEVER  IT  SAVES  A  CONSIDERA VL3  AMOUNT  CF  COMPOETER  TIME  C 
C  DURING  SIMULATION  OF  INFILTRATION  WilEM  CT  IS  SMALLEST  C 
C  IF  THIS  FEATURE  IS  NOT  DESIRABLE  THE  USER  MAY  IGNORE  C 
C  IT  uY  REPLACING  N  =  30  BY  N =CL/DZ ♦ 0 . 0 1  IN  THE  MAIN  C 
C  PROGRAM  AND  DELETING  ALL  CALL  INDTDZ  STATEMENTS  C 


^*** ************************************************************** *************c 

C 

SUBROUTINE  INDTDZ 
CCMMON/LJ/  N,NM1,NM2,NF1,NF2 
C0MMCN/L4/  ALPnA, BETA, CT.DZ 
COMMCN/L5/  NX,NX1,NRMAX,CON1 
COKMCN/LSy  TIME,TINF,TCYC 
CCMMCN/Ll  1/  CL,CL1,CLz,L1,L2 
N  =  N  ♦  1 C 

IF  (TIME.  JF..TINF)  N  =  CL/DZ*0.10 

N  M 1  =  N-  1 

KM2=NM2 

NF1=N»1 

NPc=N*2 

N  X  =  N  X ♦  10 

iF  (tlX.GT.  Nh.MAX)  N X  =  N R M A X 
NA1=NX-1 
r  e'upx 

FNJ 
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********** *************** *•*•*****•***•**••****•*****•*•••**• *************** **c 
HATER  F  L  C  H  PROGRAM 
********** ******* ******* *******************************  ****************** *****£ 


THIS  SUBROUTINE  PROVIDES  THE  SOLUTION  FOR  THE  HATER  FLON  EQUATION 
FOR  A  HOMOGENEOUS  OR 
FOR  A  LAYERED  SOIL  PROFILE 

THE  METHOD  OF  SOLUTION  IS  BASED  ON  FINITE  DIFFERENCE  APPROXIMATION 
(  SEE  TEXT  ) 

THE  BOTTOM  BOUNDARY  CONDITION  IS  INCORPORATED  IN  THE  SOLUTION 
AND  IS  APPLICABLE  FOR  A  HATER  TABLE  BOUNEARY  Oh  A  SEMI  -  INFINITE 
(I.E.  DEEP)  SCIL  PROFILE 

THE  SURFACE  BOUNDARY  CONDITION  (  FLUX  BOUNDARY  CONDITION  )  IS 
PROVIDED  FROM  MAIN  PROGRAM. 


SUBROUTINE  HAIEE 
CCMMCN/L1/  C  (310) ,Y  (310) 

COMMCN/L2/  A A  (310)  , HU  (310) ,CC (310) ,R (310)  ,RDI5T (310) 

COMMG.N/L3/  N ,  N M 1 ,  N  M2 ,  N 1 1 ,  NP2 

COMMON/L4/  ALPHA, BETA, ET,DZ 

COMMCN/L5/  NX.NX1, NRMAX,CCN1 

C0MMCS/L6/  ELNH4,PLN03,DNnRF 

COMMO N/L 7/  SFLUX,  ST  ,'.M,  QK,CSNH4  ,  CSNC3,  DISP,  XL 
COMMON/Ld/  XXX  (30)  ,C1  (3C)  ,C2  (30) ,C3  (30) ,C4  (30) 
COMMON/L9/  TIME,TI NF ,TCYC 

C0 1 MO N/L 10/  H  (310)  ,CON  (  3  1  0) , CAP  (3 1 0)  , TH  (  3 1 0) 

CCMMCN/L 1 1/  CL,CL1,CL2.L1,L2 

COM MO  N/L  1  2/  AC  1  ,BC  1  , AT  1 , ET 1 

COM MO  N/L 1 3/  Ao2,UC2,AT2,BT2 

COMMON/L 14/  AC3,bC3,AT3,ET3 

COMMON/ LI  5/  ROU1  ,THS1, ROU2 , THS2 , SOU  3, 7HS 3 

CO  1  MO  N/L  1 6/  FEXl,RNITl,RDNIT1 

CCMMCN/L 1 7/  HEX2,RKIT2, HDNIT2 

COM  MC  N/L Id/  R2X3,&NIT3,RDNIT3 

CALL  H  PROP 

DC  1  1=1, NM1 

A A  (I) =CA?  <!♦  1) ♦ ALPHA* (COS  (1*1) *CON (I) ) 

DU  (I) =-ALPUA*CGN (1*1) 

CC  (I) =-ALfriA  *CC  N  (1*1) 

1  CONTINUE 

DC  2  1=1, Nil 
X  1  =  C At  (1+1)  *K  (1+1) 

X2  =  ALPMA*CCN  (l)*ii(I)-ALPnA  +  «(I  +  1)*(CON(l+1)  *CON  (I)  ) 

X  3  =  A  L  E  il  A  *CU  N  (1  +  1) *b  (1+2) 

X4  =  -RETA*  (CON  (1  +  1) -CON (I) ) 

3  (I) =  X1  +  X2  +  X.1  +  X4 

2  CONTINUE 

r .  (1)  =h  (1)  ♦ALPHA*CON  (1)  *!)  (1) 

:  (.(Ml)  =  R  (NMl)-rL  (NM1)  *3  (NPl) 

DC  r.  I  =  1 , N  X  1 

Xrj  =  -DT*2T«CCN  C*1)*rDIST(i*1)/XL 
»  (I)  =C  (1)  *T.i  (1)/  (Th  (I)  +  XS) 
v  (1)  =Y  (I)  *Tn  m  /  (TH  (i.)  +  X9) 
nj  =  5  U)  ♦*- 
CC. .TINGS 

BALL  TMDT(AA,B2,CC,R,NM1) 

DC  !  K  =  2  ,  N 
‘  (  v)  =l«  (F-1) 

Eh!  I ;  N 
END 


non 


non 


SCIl-HATER  FhOPERTlES  PROGRAM 


c  c 

C  THIS  SUBROUTINE  PROVIDES  THE  SOIL-HATER  PROPERTIES  FOR  EACH  SOIL  C 
C  LAYER  IN  THE  SC1L  PROFILE  NAMELY  THE  HYDRAULIC  CONDUCTIVITY  AS  C 
C  A  FUNCTION  OF  HATER  CONTENT  OR  SUCTION  AND  THE  HATER  CONTENT  -  SUCTION  C 
C  RELATIONS  HIP.  FROM  THE  LATTER,  THE  HATER  CAPACITY  TERM  IS  CALOLATED  C 
C  IN  Td IS  EXAMPLE  MATHEMATICAL  EXPRESSIONS  ARE  USED  TC  DESCRIUE  THESE  C 
C  SOIL  HATER  PROPERTIES  FOP  INDIVIDUAL  SOIL  LAYERS  (SEE  TEXT  ).  C 


C 
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SUBROUTINE  *  p  ROP 

COMMCN/L2/  A  A (410) ,Bd(310) ,CC (310) ,R (JlO) ,RDIST (31 C) 
COMKCN/L3/  N,NM1,NH2,NP1,NP2 
NX,NXl,NAMAX,CON1 
H (310) ,CCN (311) ,CAP(310) ,TH  (310) 
CL,CL1,CL2,L1,l2 

ACI, EC1,AT1, BT 1 
AC2,3C2,AT2,B?2 

ACJ, BC3, AT3, BT3 
FOU1,THS1,tOU2,THS2,GCUJ,THS3 


COINON/LS/ 

COMYCN/L  10/ 

CCMM0N/L1 1/ 

COMMO  N/L 1 2/ 

COMMCN/L 1 3/ 

CCIMCN/L 14/ 

COM MO N/L IS/ 

DC  9  C  1  =  1, Li 
R  <I)=H(I) 

IF  (R  (I) .GT.O.G)  R  (I) =0-0 
XX=THS1*  (1.0*  (-B (I) /ATI) ) ** (-2) 

CAP  (I) =X X/AT 1 

TH  (I) =T  US  1/  <1.0  (-B  (I) /ATI) ) 

CON  (I)=AC1*LXP(BC1*  (id  (I) *TH  (I+1))/2) 
CONTINUE 
I  =  L  1 


CON  (I)  =AC1*LXP  ( b  C 1  *  T  H  (I)  ) 

LI  =  L  1+2 

DC  92  I =LI , L2 
E  (ij  =U  (I) 

IF (R (I) .ST. 0.0)  R  (I) =0 . 0 

XX  =  T  332*  (1.0  ♦  (-R  (1)/A1^)  )**  i-2) 

CAP  (I) =  XX/AT 2 

Til  (I)  =THS2/  (1 .0»(-R  (I) /AT 2)  ) 

CON  (I)  =AC2*EXP  (BC2*  (TH  (I)  ♦TH  (I*  1)  )  /2) 
CONTINUE 


I  =  L2 

CON  (I)  =  A  C  2  *  F  X  £'  (BC2*TH  (I)  ) 

C 

I  =  NP  1 
R  (I)  -H  (I) 

IF  (R  (I)  .  ST.O.  0)  R  (I)  =0.  C 
X  X  =  T  i!  S  3  *  (I.V.*  (-(•’  ( I )  /  A  T  3 )  )  **  (-2) 
CAP  (1) =  X X/ AT  3 

TH  (I)  *TI»3J/{  1.C»  (-;  (1) /AT  J)  ) 

CON (I) = AC 3  X  F  (BC3*T d  (I) ) 

c 

LI. =  L2»2 
DO  43  I=LJI,N 
R  (I)  =H  U) 
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IF  (R  (I) .GT.0.0)  B  (I) =0. 0 
TH  (I)  =THS3/  ( 1 .  C*  (-B  (I)  /AT3) ) 

CAP  (X)  =XX/AT 3 

XX=THS3*  (1.0*  (-B  (I)/AT3) ) •  *  (-2) 

COM (I) =AC3*£XP  (BC3*  (TH(I) *TH  (1*1) )/2) 
93  CONTINUE 
1=11*1 

TH  (I)  =  (TH  (1-1)  *TH  (I*1))/2 
COM  (I)  =  (COHO-1'  *COM(I*1))/2 
CAP(1)*(CAP IT  CAP  (1*1)  ) /2 
I=L2*  1 

TH(I)  =  (TH(l-1)*TH(I*1))/2 
CON(l)  =  (CGM(I-1)  *CON(I*1))/2 
CAP (I) *  (CAP (1-1) *CAP  (1*1) )/2 

COM 1=AC 1 *EXP (BC1*TH (1)) 

RETURN 

END 


NITRIFICATION  RATE  PROGRAM 

********* ***************************** ******* ******************* ************ 

THIS  FUNCTION  SUBPROGRAM  PROVIDES  THE  BATE  COEFFICIENT  FOB 
NITRIFICATION  AS  A  FUNCTION  OF  SCI!  HATER  SUCTION  FOB 
INDIFIDUAL  SOIL  LAYERS  (SEE  TEXT  ).  THE  USER  MAY  INCORPORATE  OTHER 
VARIABLES,  I.  E.  TIME.  SOIL  DEPTH,  ETC.,  AS  DESIREE. 

************* ******•••**••«***•*** **************** ••***••**••****•*•****•• 

E  UNCTION  Eli  (W  ,  tlh,  HC ) 

CCIflON/tV  ALPHA, LETA, Cl, DZ 
COMNCN/L9/  TIME,!  INF, TCYC 
COMMGN/L 1 1/  CL,CL1,CL2,L1,L2 
CCMMCN/L16/  F<EX1,BNIT1,RCNIT1 
COJ1MON/L17/  BEX2,KNIT2,SDNIT2 
CO’IMON/Lla/  REX3,hNI73, RDNIT3 
ZZ1=0.0 
NH=-HH 

IF  (H  I. ST.  15000)  RETURN 
Z=N*DZ 

IF (Z.GT.CL2)  GO  TO  10 
IE  (Z .07. CL  1 )  GO  TO  5 
RK 1  =  RNIT  1 
GO  TO  IS 
5  CONTINUE 
RK1=RNIT2 
GO  TO  IS 
10  CCNIINUE 
8K1=5NIT3 
15  CONTINUE 

IF  (HU. SI  .  10.0)  GO  TO  1 
RETURN 

1  IF  (Hu. GT. 50.0)  GO  TO  2 
ZZ1=8K1*  (Wil-IC.O)  *0 .005C 
RETURN 

2  IF  (NH.GT.  100.0)  GO  -C  3 

7.Z1  =  KK1*  (0.200*  (WU-50.0)  *0.0060) 

SETIIRN 

3  IF  (WH.GT. 433.0)  GC  TO  4 

ZZ1=KK1*  (O.5C0* (H3-100.C) *0.00150) 

RETURN 

4  ZZ1  =  PK1*(1.00C-(w:i-4J3.0)*O.JC020) 

IF  (Mil. UT.  1300.0)  Z ’  1  -  F  3  1  *  (0.  950-0. 050* HH/ 1000. 0) 

RETURN 

tNU 
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DENITRIFICATION  RATE  PROGRAM 

C* ****************************************************** *••***•**« ******* 

c 

C  THIS  FUNCTION  SUBPROGRAM  PROVIDES  THE  RATE  COEFFICIENT  FOB 

C  DENITRIFICATION  AS  A  FUNCTION  OF  SOIL  HATER  CONTENT  FOR 

C  1ND1FIDUAL  SOIL  LAYERS  (SEE  TEXT  ).  THE  USEE  MAI  INCORPORATE  OTHER 

C  VARIABLES,  I.  E.  TIME,  SOIL  DEPTH,  ETC.,  AS  DESIRED. 

C******** ************************************* *********************************0 

C 

FUNCTION  Z22  (M, UH, WC) 

COMMCN/L4/  ALPHA, BETA, DT,0Z 

COMHON/L9/  1IME,T1NF,ICIC 

COMHON/L  1  1/  CL,CL1,CL2,L1,L2 

CCMMCN/L 1 5/  BCU1,TUS1,BOU2,THS2,ROU3,THSJ 

CCMMON/L16/  REX1,RNIT1,RDNIT1 

CCMMON/L 17/  HEX2,6NIT2,BDN1T2 

CCMMON/Ll 8/  REX3,RNTT3, 8DNIT3 

ZZ2  =  0 . 0 

2  =  M  *  D  2 

IF  (Z . GT.CL2)  GC  TO  10 
IF  (Z.GT.CL1)  GO  TO  5 
RK2=REIN 1 
KSAT  =  TiiS  1 
GO  TO  IS 
5  6S2-RDIN2 
NSAT=1HS2 
GO  TC  15 
10  RK2  =  HD1N  3 
HSAT  =  T  HS3 
15  CONTINUE 

IF <  (WC/WSAT)  .LT.0.80)  RETURN 

z,Z2  =  RK2*  (HC-0.  80*WSAT)  /  (0.  10*RSAT) 

IF  (HC.GE.  (0. S0*NSAT) )  iZ2=RK2 

RETURN 

END 
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AMMONIUM 


E  X  C  H  A 


E  U  0  G  R  A  (1 


C  THIS  FUNCTION  SUBPROGRAM  PROVIDES  THE  RETARDATION  FACTOR 

C  FOR  AMMONIUM  EXCHANGE,  where  r  is  a  function  of  soil  hater  content  and 
C  BULK  DENSITY  OF  INDIVIDUAL  SOIL  LAYERS. 

C* ************************ ********  ****** ************** **♦♦•**** ******* *********c 


FUNCTION  ZZ3  (N, HH, HC) 

COMMON/L4/  ALPHA, BETA, DT,DZ 

COMMON/L9/  TIME, TINE, TCYC 

COMMON/L 1 1/  CL,C11,CL2,H,L2 

COMMON/LI  5/  RCU1,THS1,R0U2,THS2,R0U3,THS3 

CCMMGN/L16/  REX1.RNIT1,  RDNIT  1 

COMMON/L 17/  REX2,RNIT2,RDNIT2 

COilMCN/L  1 8/  REX3,RNIT3,  RDNIT3 

Z=M*CZ 

IF  (Z.GT.CL2)  GO  TO  10 
IP  (Z.GT.CL1)  GO  TO  5 
ZZ3=1.0«-REX1*BCU1/NC 
RETURN 

LZ3=1.0»REX2*EOU2/ilC 

RETURN 

ZZ3=1.0*REX3*rtCU3/HC 

RETURN 


o  o  r>  n  n 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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A8H0NIUH  TBANSPOBT  &  TBANSFCBHOIIONS 

P  R  C  G  B  A  H 


THIS  SUBROUTINE  PBOVIDES  THE  SOLUTION  TO  THE  AHBONIUR  TRANSPORT 
AND  TRANSFORMATION  EQUATION  UNDER  TRANSIENT  FLOB  CONDITIONS. 

IT  ALSO  CALCULATES  THE  A7HONIDH  UPTAKE  BT  PLANT  BOCTS.  THE 
HETHOD  OF  SOLUTION  IS  THE  FINITE  DIFFERENCE  APPROXIMATION  HETHOD. 
(  SEE  TEXT  ) . 

THE  HATE  OF  NITBIFIC ATION  ,  DENITRIFICATION  AND  THE  DISTRIBUTION 
COEFFICIENT  FC8  NH4-N  EXCHANGE  ABE  OBTAINED  FBON  FUNCTIONS 
ZZ1,  ZZ2,  AND  ZZ 3  ,  RES t ECTI VELT. 


SUBROUTINE  AMCNIA 
NH-4  PRGKAM 


COSHCN/L 1/ 

COMHCN/L2/ 

CCMMON/L3/ 

COHNCN/L V 

C0.110N/L5/ 

CCMMON/L6/ 

C01H0N/L7/ 

CO-THON/LV 

CCMHON/L 10/ 

CCMICN/L  1 1/ 

CCM10N/L15/ 

0C1MCN/L16/ 

CC1.1GN/L17/ 

CCHHCN/L1B/ 

M=1 


C  < J  10)  ,t  (317) 

AA  (310) ,  EB (310) ,CC  (3 1 0) , H  (3 1 0) , HDIST ( 3 10) 

N,Ntll,KH2,KEl,NP2 

ALPHA, BETA, 01, D2 

NX,NXl,NbHAX,CON1 

ELNU4, PLNC3, CNITHF 

SFLUX,EI,C7.CK,CSNH«,CSN03,DISP,XL 

TIME,TINF,TCTC 

H  (310) , CON  (310) .CAP  (3  10)  ,1H  (310) 

CL,CL1,CL2,L1,L2 

H071,THS1,H0U2,inS2,R0n3,THS3 

SEXI.BNITI, RDNIT1 

KEX2,RNIT2,3DNIT2 

1;Ea3,HNIT3,HCNIT3 


■FLX=-CON  (M)  *  (H  (?!♦  1) -H  (M)  )/DZ»CON  (.1) 

SSiNF=HFIX 

FF-DZ42.0 

C  (1)  =  (SSINP*FF*CSN1I4»DISP*TH  (1)  *C  (3)/2)/  (SSINF*FF*DISP*TH  (1) /2) 
IF  (SFLUX.LF-.C.O)  GO  TC  13 
M  =  2 

V?P=-CON  (M)  *  (H  (H*1)-:I(M)  ) /DZ»CON  (1) 


DC  5  1=1, NM1 

fi  K  K  =  7  7  3  (H  ,  K  (1)  ,TH  (H)  ) 

A A  (i) =&KK*2. 0*ALFHA*DlSP-fc21A*V?P/TH  (M) 

EF  (i) =LFTA*VIl/Trf  (tl) -ALFHA*DiSP 

B  (I) =HKK*C  (h) ♦ALPHA*31SP* (C (1*1)-2.3*C (T) *C (M-1) ) 

RK1=ZZ1  ( 7  ,  H  (  H)  ,  r:i  (tl)  ) 

E  (I)  -b  (I)  -DT*r>Kl*C  (.7) 
n=i+2 

VPP*-CGN  (H)  *  (tl  (T*  1)  -ii  (H)  )/o:*CON  (7) 

CC  (I) =-ALPHA*CISF 
CCMT  IN'JF. 

.1=1 

KK<  =  Z73(K,!I(F)  ,TU(fl)) 

A  A  (Nil)  =  EKK»ALcllA*DISP 
B (1)=B (1)  ♦ALiHA*DISF*C(1) 
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X 1-0.0 

DO  6  1*1, 8X1 

PB=Q8*BDIST  (1*1)  *C  (1*1)  /  (QK*C  (1*1)  *1  (1*1)  ) 
X1«X1*DT*P8 

6  B(I)*B(I)-Dt*PB/TH(I*1) 

PLBH4*PL8B4*X1*DZ 
60  TO  14 
C 

13  COBTIBUE 
C(1)=C(2) 

8*2 

TPP*-COB  (8) *  (H(8*1)-B(8))/DZ*COB  (8) 

DO  11  1*1,881 
BKK=ZZ3  (8,8  (8) ,TB  (8) ) 

88 (I) *BKK*2. 0* ALPHA* DISP-BET A*VPP/TH  (8) 

BB (1) *BET8*7PP/TH (8)-ALP8A*DISP 

B (X) *BKK*C (8)  *ALPH8*DISP* (C (8*1)-2. 0*C (8) *C  (8-1) ) 
BK1*ZZ1  (8,B(B),TB  (8)  ) 

B  (I) =B  (I) -DT*BK1*C (8) 

8*1*2 

TPP=-COB  (8) * (8 (8*1) -8 (8) ) /DZ*COB (8) 

CC(I) *-8LPHA*0ISP 
11  COMTX8D8 
8*8 

BKK=ZZ3 (8,8(8) ,T8(8)) 

88  (B8 1) =BKK* ALEHA*DXSP 

8*2 

TPP*-C08  (8) * (8 (8*1) -8 (8) ) /DZ*C08  (8) 

BKK*ZZ3 (8,8 (8) ,TH  (8) ) 

88 (1) =BKK*ALPH8*DXSP-BETA*VPP/TH  (8) 

C 

X1=0.0 

DO  7  1*1, 8X1 

PB=Q8*BDIST  (1*1)  *C(I*1)/  (QK*C  (1*1)  *X  (1*1) ) 
X1»X1*DT*PB 

7  B  (I)  =B  (I)  -DT*PB/TH  (1*1) 

PL8H4=PLBH4*X1*DZ 

14  C08TI80E 

CALL  TBIDB(8A,BB,CC,B,881) 

DO  IS  1*2,8 

15  C(X)*B(I-1) 

C  (BP  1)  *C  (8) 

RETOBB 

EBD 
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N  I  T  B  A  T  E  TRANSPORT  6  TBANSPORHOTIOS 

PROGRAM 


THIS  SUBROUTINE  PBOVIDES  THE  SOLUTION  TO  THE  HITBATE  TBAHSPOBT 
AMD  TBANSFOBHATIOM  EQUATION  UNDEB  TBAMSIEMT  PLOW  CONDITIOBS. 

IT  ALSO  CALCULATES  THE  MITBATE  UPTAKE  BI  PLANT  BOOTS.  THE 
METHOD  OP  SOLUTION  IS  THE  FINITE  DIFFERENCE  APPBOXIBATION  METHOD. 
(  SEE  TEXT  ). 

THE  BATE  OP  MITBIFIC ATION  ,  OENITBIPICATIOM  AND  THE  DISTRIBUTION 
COEFFICIENT  FOB  NH4-N  EXCHANGE  ABE  OBTAINED  FBON  FUNCTIONS 
ZZ1,  ZZ2,  AND  ZZ3  ,  RESPECTIVELY. 


SUBROUTINE  NITBAT 
C  N03-  PBOGBAH 

CONHON/L 1/  C (310) ,1  (310) 

C0HH0N/L2/  AA (310) ,BB{310) ,CC  (310) ,R  (310) ,RDIST(310) 

COHMCN/L3/  N,NN1,NH2,NE1,NP2 

COMHON/L4/  ALPHA, BETA, DT,DZ 

COHHOM/LS/  NX,NXl,NBHAX,CON1 

COMHON/L6/  EL NH4,PLN03,DNITBF 

COMHON/L7/  SFLUX,ET,QH, QK,CSNH4,CSN03, DISP, XL 
C0HM0N/L9/  TIHE,TINF,TCYC 

COMHON/L 1 0/  H  (310) ,CON (310) ,CAP  (310) , TH  (310) 

COMMON/Ll 1/  CL,CL1,CL2,L1,L2 

COMHO N/L 1 5/  BOU1,THSl,ROU2,THS2,BOU3,THS3 

COHMON/L 16/  REX1,RNIT1,RDNIT1 

COHHON/L18/  BEX3,BNIT3,BDN1T3 

H=1 

NFLX=-CON  (M) *  (H  (H* 1 ) -H  (M) ) /DZ+CON  (H) 

SSINF-NFLX 
PF=DZ*2. 0 

Y  (1) =  (SSINF*FF*CSH03*DISP*TH (1) *Y (3)/2)/  (SSINF*FF*DISP*TH (1)  /2) 
IF  (SFLUX.LE.0.0)  GO  TO  13 

N=2 

VPP=-CON  (M) • (H (H* 1) — H (M) ) /DZ*CON  (M) 


DO  5  1*1, NH1 

AA (I) *1.0*2. 0* ALPHA*DISP-BET  A*VPP/TH  (H> 

BB  (I) =BETA*VEP/TH (N) - ALPH A*D1SP 
B (I) =Y(H) ♦ ALPHA* DISP* (Y (H*1) -2. 0*Y  (H) *Y  (H-1)) 
BK1*Z21 (M, H (M) , TH  (M) ) 

RK2*ZZ2 (M , H (M) , TH  (M) ) 

B (I)  *B (I) *DT*BK1*C(M)-DT*BK2*Y  (N) 
X2=X2*DI*BK2*TH  (I) *Y  (I) 

M*I*2 

VPP*-CON  (H) *  (H  (M*1)-H  (M) ) /DZ*CON  (M) 

CC  (I)=-ALPHA*CISP 
5  CONTINUE 

DNITRP*DNITBF*X2*DZ 
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h*n 

AA  (N81)=1.0*ALEBA*DISP 
I  (1) =  B  (1) *ALPBA*DISP*Y  (1) 

C 

11*0.0 

DO  6  1=  1, MX  1 

PN=Q8*BDIST  (1*1)  *Y  (1*1)/  (QK*C  (1*1)  *Y  (X«- 1|  ) 
X 1=X 1  *DT*FM 

6  B  (I) =B  (I) -DT*FN/TH  (1*1) 

PIN03=PLN03*X1*DZ 
GO  TO  14 
C 

13  CONTINUE 
Y(1)*H2) 

H=2 

VPP=— CON  (H) •  (a  (8*  1)  -H  (H) ) /DZ*CON (H) 


DO  11  1  =  1,  Ntll 

Aft  (I)  =  1.0*2. 0*ALFHA*DISP-BETA*VPP/TH  (8) 

BB  (I) =BE7A*VPP/TH  (8) -ALPUA*DISP 
B (I)  =Y (8) *ALPHA*DISP* (Y (8*1) -2. 0*Y (8) *Y  (8-1) ) 
BK 1  =  ZZ 1 (B,H (8) , TH  (8) ) 

HK2=ZZ2 (B,H (8) ,TH  (8) ) 

B (I) =B  (I) *DT*BK 1*C (8) -DT*BK2*Y  (8) 
I2=X2*DT*BK2*TH (I) *Y  (I) 

8=1*2 

f PP=-CON (8) * (H (8* 1) -H (8) ) /DZ*CON  (8) 

CC  (I) =-ALPHA*CISP 
11  CONTINUE 

DNITBP=DNITBF*X2*DZ 

H=N 

AA (NB1) =1.0*ALPB4*DISP 

8=2 

VPP*-CON  (8) •  (H  (8* 1) -B  (8) ) /DZ*CON  (8) 

AA  ( 1 )  =  1 . 0 *AL  PHA*DISP-BET A*VPP/TH  (8) 

X1=0.0 

DO  7  1=1, NX1 

PN=Q8*BDIST  (1*1)  *Y  (I*1)/(QK*C  (1*1)  *Y  (I*  1) ) 
X1=X1*DT*EN 

7  B(I)=B(I)-DT*EN/TH(I*1) 

PIN03=PLN03*I1*DZ 
14  CONTINUE 

CALL  TBID8(AA,BB,CC,B,N81) 

DO  IS  1=2, N 
15  I  (I)  =B  (I-  1) 

Y  (NP1)  =Y  (N) 

BETDBN 

END 
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P  S  0  G  8  A  H  OUTPUT 


THE  FUNCTION  OF  THIS  PBOGHAH  ID  TO  PRINT  THE  BESULTS  FBOH 
HODEL  PB EDICT ION  IT  SPECIFIED  TIRES.  IN  ADDITION  IT  CALCULATES 
THE  TOTAL  AHOUNTS  OF  NITROGEN  IN  THE  SOIL  SISTER  AND  THAT  TAKEN 
OP  BY  THE  PLANTS 


SUBROUTINE  OUTPUT 
COHHON/L1/  C(310)  ,Y  (310) 

CORHON/L2/  AA (310) , BB (3 10) ,CC  (310) ,B  (310) ,BDIST  (310) 

COHHON/L3/  N,NH1,Nfl2,NE1,NP2 
COHHON/L4/  ALPHA, BETA, DT,DZ 
COHHON/LS/  NX,NX1,NBHAX,C0N1 
C0HB0N/L6/  PLNH4,PLN03,DNITBF 

C0HH0N/L7/  SFLUX, ET,QH, QK.CSNH4, CSN03, DISP, XL 
COHHON/L8/  XXX (30) ,C1 (30) ,C2  (30) ,C3  (30) ,C4(30) 

COHRON/L9/  TIRE, TINE, TCTC 

CORHON/L10/  H (310)  , CON (310) , CAP (3 10) ,TH  (310) 

CORRON/L 1 1/  CL,CL1 ,CL2, L 1 , L2 
CORHON/L 15/  BOU  1 , THS 1 , R002 , THS2 , BOU 3 , THS 3 
CORBCN/L 16/  BEX1,RNIT1,BCNIT1 
COBRON/L17/  BEX2,BNIT2,RDNIT2 
COBBON/L  18/  BEX3,BNIT3, BDNIT3 
100  FOBB AT  ( ' 1 1 ) 

200  FOBHAT  (• 1 ' ,4  0X, 'TIRE  ,  DAIS  =',F10.2/) 

299  FOBRAT  (T10, 'SOIL  DEPTH', 

$T25 , ' PBESSUBE  HEAD', 

JT43, 'SOIL-NATES  CONTENT', 

♦T65, 'HATER  FLOS'  , 

*T80, • ARHONIUR  CONCENTBATION' , 

$T 1 07, • NITBAT E  CONCENTBATION') 

301  FOBHAT (T1 4, • CH • ,T3 1 , 'CH ' ,T46, 'CH** 3/CH** 3» , T66, • VELOCITY • ,T8 2, 

••IN  SOIL  SOLUTION', T109, 'IN  SOIL  SOLUTION') 

302  FOBHAT  (T68, » CH/HH* ,T82, ' HICBOGBAHS- N/HL' ,T109, • HICROGBARS— N  /HL'/) 
499  FOBHAT (T 1 0, F8.2,T24 , FI  0. 2 ,T45,F8. 2 ,T65, F8. 4,T86, F8 . 3,T 1 1 2, F8.3) 

300  FOBHAT  (25X////, 

•25X, 'TOTAL  NO— 3  NITBOGEN  IN  SOIL  SOLUTION  PHASE  ,  HICBOGBAHS  =• , 

$F 10.3//, 

•25K, 'TOTAL  Nii-4  NITBOGEN  IN  SOIL  SOLUTION  PHASE  ,  HICBOGBAHS  =', 
IF10.3//, 

•2SX, 'TOTAL  NH-4  NITBOGEN  IN  EXCHANGEABLE  PHASE  ,  HICBOGBAHS  *', 
JF 10. 3//, 

•25X, 'TOTAL  NH-4  NITROGEN  IN  THE  SOIL  PROFILE  ,  HICBOGRANS  =  ', 
SF10.3//) 

600  FOBHAT  (25X, 'TOTAL  NITBOGEN  DENITBIFIED  ,  HICBOGBAHS  *',F10.3//) 
500  FOBHAT  (/, 

•25X, 'TOTAL  NITRATE  NITBOGEN  UPTAKE,  HICROGBARS  »',F10.3//, 

•251, 'TOTAL  ARHONIUR  NITROGEN  UPTAKE,  HICROGBARS  =',F10.3//) 

8BITE  (6,  100) 

TIHH»IlHE/24.0 
WRITE  (6,200)  TIHH 
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NBITE  (6,249) 

BBITE  (6,301) 

WRITE  (6,302) 

BFLX=SFL0X 
DC  20  I=1,NP1,2 

ZZ= (1-1) *DZ 

WRITE  (6,499)  ZZ, H (I) ,TH (I) , WFLX,C  (I) ,T  (I) 
»FLX=-CON  (1)  *  (B  (1*1)  -B  (I)  )/DZ*COM  (I) 

20  CONTINUE 
H-HP1 

DO  3C  1=1, H 
AA  (I)  =TH  ( I)  *C  (I) 

30  BB  (I)  *TH  (I)*Y  (I) 

CALL  QSF  (DZ,AA,B,tl) 

XNH4=R(N) 

CALL  CSF  (DZ,BB,B,B) 

XN03  =  B  (H) 

CALL  CSF  (DZ,C,R,NP1) 

XT  =  B  (DPI) 

CALL  CSF  (DZ,C,B,11) 

X 1-R  (LI) 

CALL  CSF  (DZ,C,B,L2) 

X12=R  (L2) 

X2=X 12-X 1 
X3=XT-X 12 

ENH4=X1*ROUl*REXl'»X2WHOU2*REX2*X3*ROU3*BEX3 

XTNH4=XNd4+ENH4 

WRITE  (6,300)  XN03,XNH4,ENH4,XTNH4 
WRITE  (6,600)  DW1TRF 

WRITE  (6,500)  PLN03, PLHH4 

RETURN 

END 


****************** **********************  ******** ***************************** *q 

C 


TRIDIAGONAL  RATRIX  PROGRAM 


THIS  SUBROUTINE  PROVIDES  SOLUTION  TO  THE  TRIDIAGCNAI  HATBIX  - 
VECTOR  EQUATIONS  FOR  SUBROUTINE  WATER,  ANONIA,  AND  NITRAT 

q******** ********************************************************* *************C 

C 

SUBROUTINE  TRIDH  (A, B,C , D, N) 

DIMENSION  A(1),B(1),C(1),D(1) 

DO  1  1=2, N 

C(I-1)=C  (1-1 ) /A  (1-1) 

A(I)=A(1)-(C(1-1)*B(I-1)) 

1  CONTINUE 

DC  2  1=2, N 

D(I)=D(I)-(C(I-1)*D(I-1)) 

2  CONTINUE 
D(N)=D(N)/A(N) 

DO  3  1=  2 ,  N 

D (N* 1-1)  =  (D (M* 1-1)- (B  (N ♦ 1-1) *D(N*2-I) ) )/A  CN* 1-1) 

3  CONTINUE 
RETURN 
END 
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APPENDIX  B:  EXAMPLE  OF  INPUT  AND  OUTPUT 


INPUT  DATA 


INITIAL  DT,  Hi.  =  0.C1C0C 
INITIAL  02  ,  CH=  1.C0C0Q 


riox  <Tr  HASTE  H  AT  EH  APPLICATION,  CH/  ft  =  C. 50000 
EVAPOT RANSPIBAT1C  N  BATE,  CF/Eh  =  C.C1CCC 

NITROGEN  LPTAKE  RATE,  BTC f COR AK-N/C t  Cf  SCOT  LENCIH  p EH  HOUR 

SICEAELIS  CONSTANT  =  1.CCC00 

CONCENT  B  ATI  CM  OF  APPLIED  NK4-N  ,  EG/L 1 1  RE  =  25.C0C.3C 
CCNCENTBATICN  OF  APPLIED  NC’-N  ,  BG/111M  =  0.0 

SOLUTE  DI SPEHSI ON  ICE  E  f 1C1 E  AT ,  CH»*2/HF  =  2.R3CJ0 


TOTAL  LENGTH  OF  SCll  EFCEILE,  CE  -  15C.CCCCC 

SOIL  DEPTH  TO  THE  EIKST  SOIL  LAYER,  C.T  =  15.C0C00 

SOIL  DEPTH  TO  THE  SEC.CSC  SCll  LAYER,  CM  =  45.00000 


SOIL  NATEF  PARAMETERS  ECR  THE  IIRST  LAYER 

SOIL  NATEf  PAPAS  ETE  RS  FOR  Tht  SECCNE  LAYER 

SC1L  NATEf  PARAMETERS  PCf  THE  THIRD  LAYEF 

FIBST  LATEH  ;  fcUIK  CEHSLTY  =  1.41CCC 

SECOND  LAYER  ;  tUIK  DENSITY  =  1.SVCCC 

THIRD  LAYRP  ;  EIILK  EF  MS  IT  Y  =  1.ESCCC 


C.9603E-0S 
0. 2200E- 05 
C.210JE-05 
SAT  U5AT ION 
SATURATION' 
SATURATION 


0. 00  ICO 


0. 27t 3E*C2 
0.30705*02 
0.  3887’J*02 
0.  4  4  C  0  0 
0.42000 
0. 34000 


FIBST  LAYER: 


SECCND  LAYER 


THIRD  LAYER: 


DURATION  d  HASTE 
SCHEDULE  CF  HAST! 
NUMBER  OF  CYCLES 


E.H4-F  T.XCHANGEAL1F  CCFFMCIEMT, 
NITPIUCATIC::  RATE  CCET.,HH-1  - 
t'MTl.lf  ICATICS  I  AT  i  COM.,  fci.-l 
NH4-N  EACH  AUG  E  A  El  E  CO  2  r  E  J.CI  ENT  , 
NiTolEICATICR  FATE-  COEF  .  ,'l  K- 1  = 
EEN1TNIMCA1  ICS  t  AT  F  CCEF.,  Hh-1 
NH4-  A  EXCUANGIAFIE  CCEFI1CIEN1, 

Ni  I  Hi  FTC  AT  LC  K  EAIF  COEF.,HF-1  « 

C  ENITRT  F IC  AT  If  E  FATE  CCEF.,  HP-1  = 

*  A 1 F  r  APPlICAUCn  ,  NFS  *  10.00000 

WAlEf  AP I  liC.’T  IC  N,  i  .  F.  CYCLE  CUE  AT  ION  = 
1 


0. 2530 J 


C  M  3/  G  M  = 

0.  13000 
=  D. 01000 
CS’/GK  =  0.250  10 

0. 10000 
=  0.01000 
C.T2/CS  -  3.25000 

0. 10C00 
=  O.CIOOO 


INci.  00000 


0.1 AOJf+OA 
0.4u03E*C2 
0.  3 0 0 0 E ♦  0 2 


0. 10 J0E+01 
C. 10002+01 
3. 100UE+01 


TINE  AT  NHICti  OUTPUT  IMS  iS 


i>  E  C  U  F  51 F  E  ,  E  i» 
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AMMONIUM  NITROGEN  UPTAKE.  MICROGRAMS  >  37.S00 


A  facsimile  catalog  card  in  Library  of  Congress  MARC  format  is 
reproduced  below. 


Selim,  H.M. 

Simplified  model  for  prediction  of  nitrogen  behavior  in  land 
treatment  of  wastewater  /  by  H.M.  Selim  and  l.K.  Iskandar.  Hanover, 
N.H.:  U.S.  Cold  Regions  Research  and  Engineering  Laboratory;  Spring- 
field,  Va.:  available  from  National  Technical  Information  Service, 
1980. 

iv,  53  p.,  illus. ;  28  cm.  (  CRREL  Report  80-12.  ) 

Prepared  for  Directorate  of  Civil  Works  -  Office,  Chief  of  Engineers 
by  Corps  of  Engineers,  U.S.  Army  Cold  Regions  Research  and  Engineering 
Laboratory. 
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